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ABSTRACT 

This  report  summarizes  the  three-year  research  program 
at  Systems,  Science  and  Software  to  investigate  the  effects  of 
meso-scale  and  small-scale  interactions  on  global  climate. 

The  research  concentrated  on  two  areas,  orographic  effects  on 
the  wind  patterns  and  effects  of  radiation  transport  on  the 
climate.  Volume  I describes  the  orographic  research  and  in- 
cludes the  theory  of  momentum  transport  due  to  mountain  ranges, 
the  formulation  of  several  computer  codes  to  calculate  the  ef- 
fects for  realistic  topography  and  wind  profiles,  and  the  ap- 
plication of  these  codes  to  various  problems  and  comparison 
with  other  calculations  as  well  as  experimental  results. 

Volume  II  describes  the  radiation  transport  research 
v^hich  produced  a benchmark  code  against  which  more  simplified 
models  can  be  compared.  This  code,  ATRAD,  is  characterized 
by  high  angular  and  frequency  resolution  and  by  the  ability 
to  calculate  radiative  atmospheric  heating  rates  taking  into 
account  molecular  absorption  and  scattering  from  arbitrary 
distributions  of  aerosols  and  particulates. 
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1.  INTRODUCTION 

The  nxamerical  prediction  of  the  general  circulation  of 
the  atmosphere  predates  most  of  the  other  applications  of  high 
speed  computers  to  physical  problems.  The  codes  which  exist 
at  several  major  research  centers  have  reached  levels  of  con- 
siderable sophistication.  These  codes  are  used  to  solve  tme- 
dependent  equations  describing  atmospheric  motion  in  a three- 
dimensional  representation.  Parametric  descriptions  are  in- 
cluded to  take  into  account  the  effects  of  insolation,  turbulent 
transport,  and  moisture. 

There  are  two  aspects  of  current  computational  capabili- 
ties, however,  which  have  been  treated  rather  crudely  in  the 
past,  and  which  have  been  the  subject  of  the  present  research 
program  at  Systems,  Science  and  Software  (S  ) : orographic  ef- 

fects and  effects  of  radiative  transfer  on  the  climatology. 
Orographic  (mountainous)  effects  are  mesoscale  phenomena,  i.e., 
occurring  on  a distance  scale  of  a few  kilometers  which  is 
smaller  than  the  typical  grid  size  in  the  Global  Circulation 
Models  (GCM) . Qualitatively,  the  effect  of  mountains  is  to 
transport  horizontal  momentum  to  high  altitudes  over  the  moun- 
tains and  long  distances  down  wind.  This  perturbation  in 
the  flow  field  can  cause  moisture  to  be  advected  to  high  alti- 
tudes and  cloud  formation  is  common.  The  clouds  and  moisture, 
on  the  other  hand,  can  alter  the  climatic  conditions  over  areas 
large  compared  to  the  mountainous  source. 

For  the  application  to  short  period  forecasts,  covering 
the  time  interval  shorter  than  a few  days,  the  details  of  the 
atmospheric  heating  by  solar  insolation  are  probably  not 
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necessary.  Over  longer  periods  of  time,  however,  the  processes 
which  transfonu  the  solar  energy  into  motion  of  the  atmosphere 
are  much  more  important  and  more  sophisticated  models  are  necessary 

The  research  effort  described  in  this  report  has  led 
to  the  development  of  computer  codes  capable  of  calculating 
the  orographic  effects  and  radiative  effects  in  considerable 
detail.  has  used  its  extensive  hydrodynamics  and  radiation 

modeling  capabilities  to  produce  codes  that  include  descrip- 
tions of  most  of  the  physical  processes  which  ure  relevant. 

Through  test  calculations  the  accuracy  of  various  physical  and 
mathematical  approximations  was  determined,  allowing  simplifi- 
cation in  the  models.  Finally,  parameterized  models  have  been 
developed  for  the  RAND  Mintz-Arakawa  two-level  Global  Circula- 
tion Model  (GCM)  which  provide  increased  accuracy  in  the  momen- 
tum transport  and  radiative  heating  calculations  at  a nominal 
increase  in  computation  time. 

The  present  report  is  a comprehensive  summary  of  nearly 
four  years  of  research.  The  editors  have  attempted  to  discard 
the  many  blind  alleys  encountered  and  reported  in  the  various 
semi-annual  reports,  and  to  condense  the  material  as  much  as 
possible,  but  the  sheer  volume  of  useful  results  has  required 
publication  in  two  volumes.  Volume  I describes  the  work  on 
orographic  modeling,  and  the  radiation  transport  work  is  de- 
scribed in  Volume  II. 
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1.1  RADIATIVE  TRANSPORT  EFFECTS  ON  CLIMATE 

In  several  intensively  investigated  problems  of  atmos- 
pheric motion,  the  role  of  radiative  transfer  is  an  important 
one.  One  of  these  is  the  global  circulation  of  the  atmosphere 
in  which  short  and  long  wave  radiation  form  a major  subsystem 
of  the  strongly  interactive  radiative-fluid  dynamic  atmospheric 
system.  In  investigations  of  climate  the  geographical  and 
seasonal  dependence  of  the  heating  of  the  atmosphere  interact- 
ing with  surface  albedo  and  cloud  cover  provides  the  source  of 
energy  for  the  atmospheric  circulation.  Another  actively  in- 
vestigated field  of  fluid  dynamics  in  which  radiation  plays  a 
role  is  the  planetary  boundary  layer.  Diurnal  changes  within 
this  layer  give  rise  to  winds  and  clouds  having  profound  in- 
fluence on  Man.  To  mention  a few  of  the  phenomena  induced  by 
the  heating  of  the  atmosphere  near  the  ground,  there  are  slope 
’:/inds,  land-sea  breezes,  heat  island  effects,  and  radiation 
and  advection  fogs.  Finally,  on  the  micro-scale  the  environ- 
ment of  plants  and  animals  is  largely  determined  by  the  radia- 
tion properties  of  their  immediate  surroundings. 

An  atmospheric  radiative  transfer  computer  code,  ATRAD,  has 
been  developed,  tested,  and  applied  to  the  calibration  of  the  ra- 
diative transfer  subroutine  of  one  of  the  climate  dynamics  com- 
puter codes.  The  code  has  the  objective  of  calculating  radia- 
tive fluxes  of  short  and  long  wave  radiation  in  order  to  de- 
termine highly  accurate  atmospheric  heating  rates.  Through 
this  formulation,  which  takes  account  of  frequency  and  angular 
dependence,  scattering  from  molecules  and  aerosols,  and  general 
surface  boundary  conditions,  we  have  made  available  a radiative 
transfer  benchmark  code  serving  as  a standard  for  comparison 
with  calculations  containing  less  comprehensive  treatments. 

The  assumptions  and  physical  considerations  forming  the 
basis  for  the  ATRAD  formulation  are  given  in  Section  2.  Also, 
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given  in  that  section  are  the  numerical  approximations  in  the 
present  code.  In  Section  3 applications  of  the  code  to  a vari- 
ety of  problems  are  presented,  including  comparisons  with  other, 
less  accurate  treatments. 
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2 . ATRAD  THEORY 


2 . 1 ASSUMPTIONS 

The  fundamental  assumptions  involved  in  the  formula- 
tion of  the  ATRAD  equations  are: 


(1) 

plane-parallel  atmosphere  (rather  than 
spherical)  and  horizontal  surface. 

(2) 

horizontal  homogeneity. 

(3) 

local  thermodynamic  equilibrium,  so  that 
true  emission  is  described  by  the  Planck 
function. 

(4) 

unpolarized  radiation  field  , 

(5) 

non-refractive  air  (rays  are  not  curved) , 

(6) 

spherical  aerosol  scatterers  of  uniform 
composition , 

(7) 

point  source  sun  (no  angular  width) , and 

(8) 

multiplying  together  the  transmissions  of 
individual  molecular  constituents  yields 
the  total  transmission. 

Assumptions  (1)  and  (5)  are  only  violated  for  large 

zenith-angle 

rays,  and  such  rays  make  only  a very  small  con 

tribution  to 

vertical  fluxes. 

Assumption  (2)  probably  represents  the  most  serious 
approximation  rn  the  code,  one  whose  impact  has  been  largely 
ignored  in  the  literature.  The  effect  of  the  most  important 
form  of  horizontal  inhomogeneity,  partial  cloudiness,  is  cus- 
tomarily accounted  for  by  taking  a weighted  sum  of  clear-air 
and  totally-cloudy  fluxes. 
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F = aF  , , + (l-a)F  , 

cj-oud  clear 

where  a is  the  fractional  cloudiness.  The  only  virtue  of 
this  scheme  is  that  it  reduces  properly  in  the  limits  a = 0 
and  a = 1 and,  therefore,  cannot  err  too  greatly  tor 
0 < a < 1 . A more  reasonable  scheme,  in  the  context  of 
ATPAD's  numerical  algorithm,  would  be  to  take  the  weighted 
sum  of  the  reflection  and  transmission  matrices  of  a partially 
cloudy  zone.  This  difficulty  could  he  resolved  with  a three- 
dimensional  Monte  Carlo  code  such  as  ^hat  of  Cox  and  McKee, 
which  should  help  to  ascertain  the  rel ati'"  importance  of 
horizontal  inhomogeneity  and  whether  it  i really  significant 
for  weather  and  climate  radiation  subroutines. 

Assumption  (3),  which  requires  sufficiently  high 
pressures  that  collisions  maintain  thermodynamic  equilibrium 
in  excited  states,  is  valid  below  about  70  km. 

Assumption  (4)  has  been  examined  by  several  investi- 
gators, for  example  Howell  and  Jacobowitz,  ^ ^ who  find  that, 
especially  in  the  presence  of  aerosols,  the  flux  errors  re- 
sulting from  a polarization-independent  treatment  are  generally 
less  than  one  percent  in  the  visible  spectrum.  Averaged  over 
the  whole  spectrum,  these  errors  are  considerably  smaller. 

[21 

A Standard  Mie  theory  treatment  of  aerosol  scatter- 
ing is  made  possible  by  assumption  (6) . Uniform  composition 
IS  not  essential,  since  layered  spherical  particles  can  be 
handled  with  some  increase  in  computational  complexity  (see 
Kerker  However,  there  is  practically  no  experimental 

data  on  layered  (water-sheathed,  etc.)  aerosol  particles. 
Neither  is  sphericity  essential  since  Mie  theory  has  been  ex- 
tended to  ellipsoids  and  circular  cylinders.  However,  solid 
aerosol  particles  (with  the  exception  of  cirrus  ice  needles) 
are  no  more  accurately  representable  as  ellipsoids  or 

*S.  Cox,  Colorado  State  University,  Fort  Collins,  Colorado, 
private  communication. 
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cylinders  than  as  spheres.  And  liquid  aerosol  particles  are 

very  nearly  spherical,  although  che  larger  ones  will  tend  to 

be  slightly  flattened  by  the  air  flow  around  them.  Further- 
. [ 4 ] 

more,  as  Hodkinson  argues,  randomly  oriented  irregular 
particles  produce  the  same  forward  diffraction,  the  same 
scattering  due  to  external  reflection,  and  the  same  scatter- 
ing due  to  the  first  refraction  as  equivalent  spheres  (spheres 
having  the  same  distribution  of  projected  area) . For  moder- 
ately absorbing  irregular  particles,  this  would  lend  to  a 
complete  equivalence  of  scattering  patterns.  Measurements 
of  Hodkinson,  ^ ^ Ellison,^  ^ ^ and  Holland  ar d Gagne ^ ^ ^ 

tend  to  confirm  the  usefulness  cf  Mie  theory  (with  equivalent 
spheres)  for  predicting  the  scattering  from  irregular 
particles  into  the  forward  hemisphere,  which  contains  almost 
all  the  scattered  photons. 

The  one  aerosol  which,  by  virtue  of  its  importance 
to  the  global  radiative  energy  budget,  f ® ^ should  probably 
not  be  approximated  as  spherical  is  cirrus  cloud. 

Jacobowitz  has  computed  the  phase  function  of  randomly 
oriented  hexagonal  ice  cylinders,  typical  of  cirrus,  and  has 
found  that  as  much  energy  is  scattered  into  the  halo,  around 
22°,  as  into  the  main  diffraction  peak.  This  effect  is  en- 
tirely missed  by  an  equivalent  sphere  model. 

The  error  incurred  by  assumption  (7) , ignoring  the 
angular  width  of  the  sun,  has  been  estimated  and  found  to  be 
negligible  as  far  as  vertical  fluxes  are  concerned. 

Assumption  (8) , that  individual  transmissions  are 
multiplicative,  is  firmly  grounded  in  experiment, and 
has  virtually  no  detectable  experimental  error  associated 
with  it. 
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2.2  BASIC  RADIATIVE  TRANSI’ER  EQUATIONS 

The  radiative  transfer  problem  in  the  Earth's  atmos- 
phere reduces  to  the  solution  of  the  seemingly  simple  equa- 
tion^^^^ 


dl 

V 

ds 


J 

V 


K 


V 


I 

V 


(2.1) 


which  states  that  the  monochromatic  (frequency  v)  radiant  in- 
tensity I^ , in  traversing  the  element  of  length  ds  , will  be 
augmented  by  sources  in  the  amount  ds  and  diminished  by 

extinction  in  the  amount  I^  ds  . In  general,  the  radiant 
intensity  and  , the  source  function,  depend  on  both  a 

spatial  coordinate  z"  and  an  angular  coordinate  (direction) 

^ at  the  point  r , as  well  as  upon  the  frequency  v . The 
time  dependence  of  these  quantities  is  ignored  because  the  j 

radiative  state  of  the  atmosphere  is  established,  for  all  j 

practical  purposes,  instantaneously,  k is  the  extinction  i 

• • ^ I 

coefficient,  which  describes  the  relative  depletion  in  th«.‘  ! 

intensity  of  the  beam,  dl^/I^  , upon  traversing  the  element  ! 

of  distance  ds  . is  in  general  the  sum  of  an  absorption  j 

part  and  a scattering  part.  describes  the  additions  made  I 

to  the  beam  intensity  along  ds  by  thermal  (Flanckian)  emis-  | 

sion  and  by  scattering.  I 

In  ATFAD's  plane-parallel  geometry,  which  is  illustrated  j 

in  Figure  2.1,  r is  replaced  by  z,  which  is  measured  from  i 

the  top  of  the  atmosphere,  and  ^ is  represented  as  usual  by  | 

the  spherical  coordinates  y = cos0  and  0 . Using  Kirch-  j 

hoff's  Law  to  obtain  the  source  function  , the  radiative  I 

transfer  equation  (2.1)  becomes  in  this  geometry  I 

2-4  ; 

I 

! 

i 


I 
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\ 


////////////  (surface) 


Figure  2.1  — ATRAD  coordinate  system 


31 

“3^  = 


V 


*'4  IT 


(2.2) 


where 


+ 6^  = volume  extinction  coefficient, 

= volume  absorption  coefficient,  , corrected 
for  stimulated  emission  (a  correction  usually 
ignored  in  the  atmospheric  radiation  literature) , 

= . ,-hvATj  _ 

T = temperature, 

3^  = volume  scattering  coefficient, 

= Planck  function 

2hv  Vc^ 
hv/kT 

G "1 
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and  is  the  phase  function,  defined  so  that 


4tt 


is  the  probability  that  a photon  entering  a voliirae  element 
around  z from  direction  will  be  scattered  into  the 

cone  df^  of  directions  around  . Since  absorption  is  ex- 
plicitly represented  in  Eq.  (2.2),  the  above  probabilities 
must  sum  to  one 


P 


V 


4tt 


1 


(2.3) 


In  what  follows,  it  is  convenient  to  break  the  absorp- 
tion coefficient  into  a "line"  part  and  a "continuum"  part, 


a. 


line  , 

% + 


cont 
a . 


(2.4) 


varies  "slowly"  with  frequency  and  includes,  for  example, 
the  water  vapor  8-12y  continuum  and  the  aerosol  absorption; 

^line  rapidly-varying  part,  due  to  absorption  lines  of 

atmospheric  constit' lents . Since  a continu’om  is  often  due  to 
the  wings  of  distant  lines,  the  separation  is  to  a certain  ex- 
tent arbitrary.  Therefore,  we  shall  adopt  the  convention 
that,  with  the  exception  of  aerosol  absorption,  continua  exist 
only  in  the  gaps  between  absorption  bands. 

Because  a volume  element  in  the  atmosphere  exhibits  the 
same  scattering  pattern  no  matter  what  the  direction  of  incidence 
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i 


, the  phase  function  depends  only  on  the  scattering  angj.e 


P 


P 

V 


(z,Ms) 


where 

y = COS0  = yy  ' + /l-y ^ /l-y  ' ^ cos  ((};-(})') 

In  order  to  azimuthally-average  the  radiation  transport  equa- 
tion (see  below),  it  is  necessary  to  azimuthally-average  P^ : j 

P^(z.U,li')  = ^ f P (z,w Jd*  i 

*^c  i 

i 

r 
( 

P (z,yy'  + /l-y^  A-y ' ^ cos(^))d9  . (2.5) 

^ I 

f 
i 

! 

The  second  step  follows  from  the  periodicity  of  the  cosine.  \ 

As  a result,  P^  does  not  depend  on  (f)'  . J 

i 

For  the  computation  of  vertical^  fluxes  and  radiative  | 

heating  rates,  the  retention  of  the  (fi-dependence  of  is 

unnecessary.  It  can  be  eliminated  by  (})-averaging  Eq.  (2.2) 

^2tt 

(operating  on  both  sides  with  J dp)  to  obtain  | 
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U3—  = o'  B + 
a Z V V 


I^(z,u')du'  - K I (2.6) 
^ V V V 


where 


(z,h)  = 


= ^ f 

Jo 


Knowing  , the  vertical  flux  F , is  obtained  from 


(2)  = j p I (z,i^)dfl  = 2ti 

J4v  ^ J_^ 


y I.^(z,y)dy 


(2.7) 


From  F^  the  total  (spectrally-integrated)  flux 
computed: 


F can  be 


F(z)  = 


= I 


V . 

min 


Fy (z) dv 


where  ^\in' '^max^  includes  all  spectral  intervals  in  which 
F^  IS  non-negligible.  From  F , the  radiative  heating  rate 
at  level  z can  be  calculated  from 


= l_  ^ 

^ pCp  dz 


(2.8) 


P -•  p(z)  = air  density 

- specific  heat  of  air  at  constant  pressure 
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Using  the  hydrostatic  approximation-  this  becomes 


3T  _ g dF 

" C"  dp  • 

Rather  than  compute  point-values  of  heating,  ATRAD  computes 
the  average  heating  of  a zone  in  terms  of  the  bound- 

ary fluxes: 


dp 


_ g 

~ ■ S 


Pi+1  - Pi 


(2.9) 


(the  slight  vertical  variation  of  g has  been  ignored) 

The  radiative  transfer  equation,  Eq.  (2.6),  is  a 
phenomenological  equation  requiring  input  data  of  three  dif- 
fer<2nt  types: 


(a)  atmospheric  structure  data 

(b)  boundary  condition  data 

(c)  absorption  data  for  relevant  atmospheric 
constituents . 

i 

The  next  three  sections  describe  ATRAD 's  needs  in  each  of  these 
areas . 
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2.3  SPECIFICATIONS  OF  ATMOSPHERIC  STRUCTURE 

The  atmospheric  structure  data  needed  to  solve  Eq. 

1 2. 6)  , with  the  use  of  Mie  scattering  theory  for  aerosols, 
consists  of  the  pressure  p,  the  temperature  T,  the  water 
vapor  density  the  ozone  density  and  the  aerosol 

(including  cloud?  number  density,  size  distribution,  and  compo- 
sition. These  must  all  be  given  as  a function  of  height. 

The  atmospheric  absorbers  other  than  water  vapor  and  ozone 
are  assumed  to  be  uniformly  mixed. 

The  details  as  to  how  the  structure  data  are  sup- 
plied to  ATRAD  are  not  given  here.  A variety  of  analytic 
aerosol  size  distributions  are  available  as  options,  or  the 
user  may  supply  his  own  size  distributions  in  card  form. 
Similarly,  a variety  of  aerosol  refractive  indices  are  avail- 
able as  options,  e.g.  water,  ice,  sea-salt,  dust,  etc. 

2.4  BOUNDARY  CONDITIONS 

The  boundary  condition  at  the  top  of  the  atmosphere 
is 


1.,  (0 , y , <J')  = S 6 (tl  - ^ ) for 

V Vo 


0 < y ^ 1 
0 < 4)  < 2tt 


where  S^  is  the  solar  flux  at  frequency  v at  the  position  of 
the  Earth  and  is  the  direction  of  the  solar  beam  in  ATRAD' s 

coordinate  system  (Figure  2.1).  Taking  the  <|)-average, 

S.. 


" Tn  " ^o^  0 < y 1 


(2.101 
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I 

I 

I 

I 

I 


I 


I 


the  cosine  of  the  solar  zenith  angle,  is  calculated  follow 
ing  the  prescription  of  Woolf from  user-supplied  values  of 
latitude,  longitude,  day  of  the  year,  and  Greenwich  time.  The 

ri5  1 

best  measurements  of  are  those  of  Thckaekara,  et.al.,^  ^ 

which  have  been  included  in  the  code.  (Actually,  the  values  o 

00 

dX  = r S dv 

J\j  ^ 

are  given  by  Thekaekara,  and  the  code  obtains  from  these  the 
values 


S.  - - — f S dv 
Av  Av  /,  V 
*'Av 

for  use  in  each  frequency  group  Av.)  is  adjusted  accord- 

ing to  the  earth-sun  distance,  v/hich  is  calculated  from  the 
day  of  the  year. 

The  boundary  condition  at  the  surface  requires , for 

its  complete  specification,  the  ground  temperature  T and 

the  bidirectional  reflectivity  ^ beam  of  intensity 

^v,inc  incident  on  a surface  from  direction  and  the 

resultant  reflected  intensity  at  angle  is  I ^ (h  . ) 

r 16l  ^ v,ref  r x ' 

then  p " is  defined  such  that^  ^ (see  Figure  2.2). 


(2.11) 


By  this  definition,  is  of  finite  order,  barring  specular 

reflection,  for  is  of  differential  order  with  respect 

^v,inc*  there  is  a specular  component  of  reflec- 

tion, we  separate  p"  into  a diffuse  part  p"  , 

V ^ ^v,d 


and  a 
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specular  part  , 


P 


II 

V 


P 


II 

v,s 


where 


:r-  P 

COS0,  V 


(2.12) 


n_  6 1 

Some  authors,  for  example  Siegel  and  Howell  , do  not 
include  the  factor  of  it  in  Eq.  (i.ll).  We  find  it  convenient 
because  v;hen  the  reflection  is  diffuse,  p^  reduces  to  an 
albedo  (flux  ratio) , the  customary  measure  of  surface  reflecti- 
vity. 


Figure  2.2.  Geometry  of  reflection. 


VJhen  radiation  is  incident  from  all  the  reflected 

intensity  along  is,  by  Eq.  (2.11), 


1 /•' 

^ (f)  ) = ± I dd>.  1 dy  . p"  (SL  ,^  ) I . (S2.  )y  • 

,ref  ' r'  n I i | i '■  i'  r'  v,inc  i 

^ 0 J Q 


Adding  on  the  thermal  emission  from  the  ground  and  phrasing  the 
whole  result  in  terms  of  the  ATRAD  coordinate  system  (Figure  2’.1), 
the  surface  boundary  condition  becomes 


I^,(Zo,y,(())  = e^(  |y  l)B^(Tg)  + ( jy  !)  (z„ , |y  | 

.2TT 


(2.13) 


-211  -1 

ij  d({.'  I dy'p"  ^(iyj,y',(()-4)')I^(z^,y',(()')y' 
0 •'0  ' 


for 


- 1 £ y < 0 
0 < 4>  < 2tt 


where  is  the  surface  level,  and  e^,  the  directional 

emissivity,  has  been  written  as  a function  of  angle  of  emis- 

r 1 61 

sion.  •'  The  absolute  value  signs  on  y in  e',  p‘  , and 

V V ^ s 

p^  ^ are  merely  because  these  surface  properties  are  customarily 
measured  or  calculated  with  the  convention  0 £ 9^'  £ 90®, 

whereas  ATRAD 's  coordinate  system  takes  0 £ < 90°  and 

90°  < 9^  1 180°. 

It  should  be  noted  that,  in  Eq.  (2.13),  (j)  is  entirely 

missing  from  the  argument  list  of  p'  and  occurs  in  the  com- 

V ^ s 

bination  cJ)-4)'  in  the  argument  list  of  p^  In  both  cases, 
this  is  justified  because  the  surface  is  isotropic,  that  is. 
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its  reflecting  properties  do  not  change  as  we  rotate  the  co- 
ordinate systein  in  Figure  2.2  about  the  vertical  axis. 

Taking  the  (^'“average  of  Eq.  (2.13)  leads  to  the  final 
form  of  the  surface  boundary  condition: 


I^(Zg,p)  = c;,  (T^)  + , ly  \) 


V " ' V g '^v,s  ' V 0 


/ 


1 


(2.14) 


+ 2|  dy'  ( ly  1 ,y  ' )I^  (Zq  ,y  ' ) y ' 

0 ' 


for  - 1 < y < 0 , where 


-2^ 

J 0 


-2TT 

‘jwf 

^ n 


(the  second  step  follows  from  the  2iT-periodicity  of  p"  ^ in 
its  third  argument) . 

There  are  three  simplifications  of  the  surface  boundary 
condition  (2.14)  which  are  used  in  almost  all  of  the  ether  ex- 
tant atmospheric  radiation  codes: 

(a)  e;  = 1,  p;^^  = p;^^  = 0.  This  is  used 

practically  without  exception  by  codes 
dealing  with  the  infrared  (IR)  spectrum, 
ignoring  the  very  real  frequency-and 
angle-variations  of  for  natural  sur- 

faces. Furthermore,  except  for  water, 
snow,  and  ice,  the  average  IR  emissivity 
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of  natural  surfaces  is  known  to  deviate 
considerably  from  unity- Nevertheless, 
the  mathematical  consequences  of  using  a 
non-unit  emissivity  (impiyinq  d ^ 

are  disastrous  for  codes  using  a simple 
quadrature  of  the  integral  transport  equa- 
tion, for  then  they  can  no  longer  frequency- 
average*.  the  boundary  term. 

(b)  p"  - = 0,  p'  given  by  Fresnel  formulae 

V / Cl  V ^ o 

for  reflection  at  a plane  interface,  v so 
large  that  B^(Tg)  is  negligible.  This  is 
used  primarily  by  codes  dealing  with  single 
frequencies  in  the  solar  spectrum  with 
water  as  the  underlying  surface. 

(c)  p'  =0,  p"  , = constant  (independent  of 

v,s  v,d 

angle  of  incidence  and  angle  of  reflection) , 

V so  large  that  B^(Tg)  is  negligible. 

This  is  used  primarily  by  codes  dealing  with 
single  frequencies  in  the  solar  spectrum 
with  a solid  underlying  surface.  p^  ^ is 
the  albedo,  and  the  reflection  is  assumed 
diffuse . 

We  believe  that  it  is  important  to  include  a fully  general 
boundary  condition,  Eq.  (2.14),  as  ATRAD  does.  Research  to 
date  suggests  that  such  factors  as  the  angular — and  frequency — 
dependence  of  e'  or  p"  , cannot  be  ignored. 

While  data  with  respect  to  these  dependences  are  limited,  they 
exist  in  sufficient  quantities  to  warrant  comparisons  with 
calculations  using  frequency  and/or  angular  averages. 
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Tv7o  theoretical  models  of  rough-surface  reflectivity 
seem  sufficiently  realistic  to  v/arrant  inclusion  in  ATRAD. 

The  first  is  that  for  a rough  v;ater  surface,  in  which  it  is 
assumed  that  the  surface  consists  of  planar  elements  each  of 
which  is  a specular  reflector.  The  slope  distribution  of  the 
planar  elements  as  a function  of  average  surface  wind  speed 
can  be  taken  from  the  experimental  work  of  Cox  and  Munk.^^"^^ 
Even  with  this  model,  the  inclusion  of  shadowing  and  multiple 
reflections  among  planar  elements  is  difficult.  Therefore 
these  effects  can  probably  only  be  included  approximately. 
Experimental  measurements  of  sea  surface  bidirectional  re- 
flectivities ^ 
tions . 


should  prove  useful  in  guiding  the  approxima- 


The  second  theoretical  model  of  rough-surface  reflect!-  j 

vitv  is  for  land  areas,  and  assumes  that  the  surface  consists  | 

of  planar  elements  each  of  which  is  a diffuse  reflector  (of  given  ' 

albedo).  This  model  is  due  to  McStravick . ^ The  slope  distri-  i 

bution  of  the  planar  elements  is  presumed  Gaussian.  As  with  the 
first  model,  shadowing  and  multiple  reflection  can  probably  be 
included  only  approximately. 

In  practice,  the  full  bidirectional  reflectivity  of  the  f 

earth's  surface  is  not  available  experimentally,  nor  is  it  j 

likely  to  be  available  in  the  neai  future,  as  a function  of  | 

latitude  and  longitude.  Therefore  the  two  theoretical  models  j 

described  above,  particularly  the  water-surface  model  applying  j 

to  so  much  of  the  earth's  surface,  assume  special  impo.-tance  | 

in  filling  this  void.  However,  ATRAD  also  makes  provis.'.on  to 
use  certain  degraded  forms  of  the  bidirectional  reflectivity 
and  directional  emissivity  which  are  more  readily  available 
fiom  experiment.  These  are  (refer  to  Figure  2.2  for  notation): 

(a)  the  directional-hemispherical  reflectivity, 
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SSS-R-75-2556 


p I ^'^2  ) = reflected  flux 

^ incident  flux  from 


r I ^ (Q  ;n.)cOSd  dfl 

v,ref  r'  1 r r 


“if  p"(h.,'ii  ) 

TT  / 1 ' r ' 


COS0  d.Q 

r r 


fl 


dp^ 


(b)  the  hemispherical  reflectivity  (albedo) , 


reflected  flux 
incident  flux 


_ •'Q 


/* I . ) COS0  . df2 . 

I v,inc  1 1 1 


(c)  the  hemispherical  emissivity. 


e . = 


emitted  flux 


black  body  flux  at 


= IQ 


7TB  (T  ) 
V ' g' 


= ? d!) 


■=/' 

A 


e^(y)y  dy 


(2.15) 


(2.16) 


(2.17) 
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The  solid-angle  integrals  in  (a)-(c)  are  over  the  upward  hemis- 
phere. Note  that,  by  the  isotropy  of  the  surface,  really 

only  depends  on  Also,  the  above  quantities  are  not  all 

independent.  By  Kirclihoff's  uav/  for  surfaces  ^ , 

1 (2.18) 

so  that,  if  either  or  is  known,  the  other  is  deter- 

mined . 

If  only  the  directional-hemispherical  reflectivity 
is  given,  nothing  can  be  said  about  the  angular  dependence  of 
the  reflected  radiation.  Therefore  ATRAD  assumes  diffuse  re- 
flection. which  means  p^  in  Eq.  (3.15)  is  independent  of 
and  is  equal  to  p^, 

p;;(yr'Pi)  = p;(Pi) 

If  rather  than  p^  is  known,  Eq.  (2.18)  is  used  to  obtain 

p^,  and  therefore  "p^. 

Note  that  the  albedo  p^,  Eq.  (2.16),  is  not  an  intrinsic 
surface  property  but  depends  on  the  incident  radiatio.i  field 

^v,inc'  ^ familiar  to  experimenters.  ^20]  when  p^  is 

independent  of  angle  does  the  albedo  become  independent  of  in- 
tensity, in  which  case  Eq.  (2.16)  reduces  to 


When  only  an  albedo  is  available,  ATRAD  uses  the  approximation 
that  p^  is  independent  of  both  and  y^  , leading  to 
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To  obtain  the  emissivity,  observe  that  if  the 
i-otropic  or  xf  is  independent  of  angle, 


emission 

then 


is 


P + 

V 


c 

V 


1 


from  Eq.  (2.18). 
he  obtained  from 
isotropic  emission 


Since  the  second  assumution 
Py.  ATRAD  then  invokes  the 
to  obtain 


holds,  can 

assumption  of 


e ' = 

V 


e 

V 


rrom  Eq.  (2.17). 

as  i.  ~ “ available,  particularly 

1.  relates  to  remote  sensing  from  aircraft  and  satellites  Much 
more  such  data  remains  classified.  However,  as  yet  no  compre- 
sive  tabulation  of  albedos,  bidirectional  reflectivities  etc 

will  be  obtained  as  needed. 


2.5 


TRANSMISSION  FUNCTIONS 


and  absorption  coefficients  have  been  calculated 

and  measured  during  the  last  two  decades  for  all  of  the  important 
tmospheric  absorbers.  However,  it  is  computationally  infeasibL 

in  individual  lines 

da.  a in  th  , ^■^®'3“'=ndy-averaged  absorption  coefficient 

da .a  in  the  form  of  transmission  functions  are  employed. 

In  terms  of  the  optical  depth  for  line  absorption 


,,ab<^i''=i>  = f ' (2)d2, 

2 1 


(2.19) 


the  transmission  function  is 
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T 


Av 


(Zj;Z2;P) 


1 /•  (Z  W22  )/y 

1 ^ V , ab 


dv 


(2.20) 


The  frequency  interval  Av  in  general  includes  many  lines. 

The  computation  of  transmission  functions  has  a long 
history.  The  earliest  attempts  v/ere  based  on  band  models,  in 
v/hich  simple  analytical  representations  of  line  strengths, 
positions,  and  shapes  were  assumed.  As  accurate  line-by-line 
absorption  data  has  become  available,  both  from  theory  and  ex- 
periment, transmission  function  computations  have  incorporated 
it.  However,  such  detailed  line-by-line  transmission  function 
computations  are  very  expensive  in  terms  of  computer  time.  Con- 
sidering that  integration  steps  as  small  as  0.001  cm  may  be 
required,  and  that  the  region  of  significant  absorption  extends 
from  10,000  cm" ‘ (It)  to  250  cm"^  (40p),  with  only  a few  gaps, 

[2  2 J 

the  magnitude  of  the  problem  is  apparent.  As  an  example,  Kyle 
used  15  minutes  of  CDC  6600  time  to  compute  transmission  func- 
tions between  a single  pair  of  atmospheric  levels  Zj  and  Zj^ 
for  a single  value  of  y,  for  the  wavelength  interval  1.7y  - 20y 
with  Av  = 1 cm~\  Multiply  this  by  the  number  of  angles  and 

by  the  number  of  pairs  of  levels  ^N^(N^-l),  and  the  computing 
time  to  obtain  a complete  set  of  transmission  functions  becomes 
truly  formidable  (for  4 angles  and  20  levels  Kyle  would  reouire 
190  hours).  While  it  is  foreseeable  that,  with  ingenuity,  such 
computing  times  could  be  reduced  substantially , it  is  apparent 
that  for  the  present  ATRAD  cannot  afford  to  calculate  transmis- 
sion functions  line-by-line  for  each  new  problem.  Therefore, 
a less  exact  method  was  sought. 

A rather  thorough  literature  survey  revealei  numerous 
^ess  exact  methods  for  computing  transmission  functions.  These 
ranged  from  simple  parameterized  band  models  such  as  Elsasser's 
and  Goody's, through  more  sophisticated  band  models  such  as 
Wyatt's, to  approaches  such  as  those  of  Smith ^ or 
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McClatchey ^ ' involving  an  empirically  determined  function  as 
well  as  empirically  determined  parameters.  The  errors  in  the 
methods  of  Smith  or  McClatchey  are  almost  within  the  uncer- 
tainties in  the  line-by-line  calculations  themselves,  and  in- 
volve orders  of  magnitude  less  calculation.  The  method  of 
McClatchey  was  selected  for  inclusion  in  ATRAD  because  of  its 
comprehensive  treatment  of  all  known  atmospheric  absorbers,  its 
accuracy,  and  its  computational  economy  (only  tabular  interpola- 
tion is  required  to  obtain  transmissions) . 

The  method  of  McClatchey,  et.al.,  as  furnished  to  us  in 
subroutine  LOTRAN  (for  LOw-resolution  TRANsmission  functions), 
contains  separate  values  of  the  transmission  for:  (a)  H2O 

vapor;  (b)  the  uniformly  mixed  gases  CO^,  N2O,  CH^,  CO,  and  O2; 
(c)  O^;  (d)  the  H2O  8-13u  continuum;  and  (e)  the  N2-N2  4p  con- 
tinuum. In  ATRAD,  the  two  continue  (d)  and  (e)  have  been  placed 
in  subroutines  separate  from  LOTRAN  because  their  transmissions 
are  simple  exponentials  (for  small  enough  Av)  and  they  can  be 
added  to  cont,  the  continuum  part  of  the  absorption  coef- 
ficient. 

LOTRAN  furnishes  transmissions  averaged  over  20  cm~^ 
intervals  centered  from  350  cm~^  to  50,000  cra~^.  To  obtain 
transmissions  over  larger  wavenumber  intervals,  it  follows 
rigorously  from  the  definition  (2.20)  that  the  20  cm"^  trans- 
missions may  be  averaged;  e.g., 

"^400-440  ""  2 ^'^400-420  "^420-440^ 

Hence  ATRAD  may  use  wavenumber  intervals  which  are  any  multiple 
of  20  cm  subject  to  restrictions  on  the  size  of  Av  . For 

wavenumber  intervals  centered  below  350  cm~^,  the  Goody  random 
band  model  for  the  H2O  rotational  band  is  used  in  LOTRAN,  on 
the  recommendation  of  one  of  McClatchey 's  colleagues  (J.E.A. 
Selby) . Parameters  for  this  band  model  were  taken  from  Goody, 
Table  5.5. 
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In  order  to  account  for  atmospheric  slant  paths  along 

which  the  pressure  p and  temperature  T,  and  therefore 
line 

“v  ' vary  that  is,  in  order  to  perform  the  integration  in 
Eq.  (2.19)  LOTPIAN  employs  a version  of  the  scaling  (or  effec- 
na£s)  approximation.  This  approximation  is  an  attempt  to 
sum  up  zi,  za  and  the  p-T  variation  of  a^^^^(z)  into  a 
single  argument  u called  the  'effective  absorber  amount.' 
u is  calculated  using  certain  a priori  assumptions  concerning 
the  pressure  and  temperature  variation  of  with  parameters 

chosen  to  best  fit  the  real  data.  Following  Goody, the  scal- 
ing approximation  assumes  the  following  decomposition  of 

line  , . , . 

(the  absorber  density  is  separated  out  because  the  mass 

absorption  coefficient  intrinsic  property,  is  preferred 

by  spectroscopists) . The  stimulated  emission  correction  has  been 
omitted  from  explicit  consideration,  but  the  effects  of  stimu- 
lated emission  are  present  in  the  LOTRAN  model  insofar  as  certain 
parameters  are  determined,  at  least  in  part,  from  measured  trans- 
missions*. The  absorption  optical  depth  (2.19)  then  becomes 

XZ2 

p^(z)aj,  (p,T)dz 
“ 1 

= u(z, ,Zj)  (2;21) 


*McClatchey,  private  communication. 
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where 


= Uj  {v)ct2  (Po/T(,  j 


u (z^ ,z^ ) 


(z) 


ct?  (P,T) 

^2  (P  0 ' 0 ^ 


dz 


(2.22) 


and  Py  and  are  taken  in  LOTRAN  as  STP  (1  atm. , 273^K) . 

u is  the  equivalent  absorber  amount  on  the  vertical  path  be- 
tween Zj  and  Zj.  The  further  assumption  is  then  made  that 
a 2 has  the  form 


(p,T)  = Ap^ ^ 


Different  values  of  Yj,  chosen  to  give  a best  fit  to  line 
by-line  transmissions,  are  used  in  LOTRAN  for  each  category 
of  absorber  (a)-(e). 


Unfortunately,  because  of  lack  of  good  data  on  line 
strengths  as  a function  of  temperature,  it  has  not  been  pos- 
sible to  obtain  best  values  of  in  the  same  way.  An 

attempt  has  been  made  to  include  the  temoerature  dependence 
of  the  line  half-width  by  using  p/T  ‘ instead  of  p,  so  that 
at  present  Y2  = ” Yj/Z,  but  it  is  not  at  all  certain  that 
this  leads  to  less  error  nor  that  T ^ is  indeed  the  correct 
temperature  dependence  of  the  half-width.  A prescription  for 

putting  temperature  dependence  into  T.,,  has  been  proposed 

[27] 

by  Kondratyev  and  Timofeev.  These  authors  claim  that 

transmission  errors  of  up  to  .1  - .2  may  result  if  tempera- 
ture dependence  is  not  accounted  for,  even  if  one  uses  the 
Curtis-Godson  approximation  instead  of  scaling.  Therefore 
more  work  is  needed  on  this  question. 
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Many  IR  radiative  transfer  studies  emplov  the 
Curtis-Godson  approximation  (see  Armstrong  ^ for  an 
excellent  discussion)  rather  than  the  scaling  approxi- 
mation. Very  few  comparisons  have  been  made,  however, 
of  Curtis-Godson  and  scaling  with  exact  calculations . 

Tv'o  such  comparisons  are  that  of  Kondratyev  and  Timofeev, 

discussed  in  the  last  paragraph,  and  that  of  Zdunkowslci 
[291 

and  Raymond.  The  first  concludes  that  Curtis- 

Godson  and  scaling  incur  comparable  errors.  The  second 
concludes  that,  while  Curtis-Godson  has  the  edge  in 
accuracy,  the  scaling  approximation  is  accurate  to  at 
least  two  decimal  places  for  a wide  range  of  path 
lengths.  Thus,  it  appears  to  be  justified  to  use  the 
scaling  approximation,  especially  for  applications  such 
as  the  present  one  in  which  errors  in  any  one  frequency 
interval  are  mitigated  by  the  integration  over  all  fre- 
quencies. 

2.5.1  Exponential-Sum  Fitting 

The  method  of  Cantor  and  Evans for  exponential- 
sum  fitting  is  now  incorporated  into  ATRAD,  although  nu- 
merous modifications  and  simplifications  have  been  intro- 
duced in  adapting  the  method  to  the  fitting  of  transmission 
functions.  The  method  is  briefly  outlined  in  the  follow- 
ing paragraphs. 
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We  desire  to  fit  a transraiasion  function  T^^(u), 
whose  values  "^'/\\)^^n^  are  gi.ven  at  a set  of  points 


0=Uj<U2<...<U  , 


with  a positively-weighted  sum  of  exponentials, 


-k . u 


T*..( 


u,  ^ , a, 


. > 0 


The  fit  is  to  be  best  in  the  least-squares  sense,  that  is, 
it  is  to  be  such  that  the  residual 


IN  t’i  -|2 

= L ”n  - L 

n=l  i=l 


is  rainimized.  The  quantities  W^  are  positive  weights  which 
may  be  assigned  arbitrarily.  In  order  to  ensure  that  small 
and  large  values  of  the  transmission  receive  equal  treatment, 
the  weights  have  been  taken  as 


SO  that  R is  in  fact  a sum  of  squared  relative  errors, 

The  u 's  are  selected  to  be  equally  spaced, 
n 


u = (n-l)Au 
n 
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Then  R ;aay  be  written 

N M . 2 

R = v;  T (u  ) - a.o”  ^ 

/ j n Av  n 1 1 

n=l  i-=l 

where 

-k . Au 

6i  = e ^ 

Let  us  presume  the  following  ordering  of  the  6^'s: 

O<e<0  <...<0w<l 

— 1 2 J>1  — 

The  problem  of  minimizing  R is  partDy  linear  (finding 
aj,...,aj^^)  and  partly  nonlinear  (finding  0 ^ , . . . , 0j^j)  . The  algo- 
rithm which  we  shall  describe  below  treats  the  linear  and  non- 
linear parts  of  the  problem  separately,  although  we  shall  see 
that  there  is  some  coupling. 

Given  a set  of  0j,...,6j^,  we  shall  call  the  procedure 
of  finding  the  aj  , . . . ,aj^j  which  minimize  R with  no  constraints 
"solving  the  linear  problem."  Since  the  a^'s  are  unconstrained, 
solving  the  linear  problem  may  produce  one  or  more  negative  a^'s. 

The  extreme  (and  well-documented)  ill-conditioning  of 
the  linear  problem  is  side-stepped,  as  far  as  the  calculation 
of  the  residual  is  concerned,  by  the  use  of  divided  differences, 
which  is  made  possible  by  the  use  of  equally-spaced  u^'s  . 

During  the  iteration  between  linear  and  non-linear  problems, 
the  a.  ■ are  of  importance  only  in  determining  which  exponential 
factors  0^  are  to  be  dropped  because  of  a negative  coefficient. 
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To  find  the  best 
that  to  the  approx ijiiant 


9.  's 
1 


we  proceed  iteratively. 


Imagine 


■1=1 


„n-l 
a . 0 . 
a.  1 


a . > 0 
1 


we  add  a new  tern  with  exponential  factor  0 and  coefficient 
a.  We  want  to  see  how  the  residual  changes  as  we  increase  a 
from  zero;  therefore,  we  evaluate  the  slope  of  R at  a=0: 


= - 2 


Z)  «nkv<%> 
n=l 


-E 

■i  =1 


1 


gn-l 


(2.23) 


which  is  simply  a polynomial  in  6.  Since  we  want  to  minimize 

the  sum  of  square.s,  we  pick  9 so  that  this  exprt?.ssion  is  as 

negative  as  possible  (if  it  is  never  negative  for  0 ^ 6 £ 1, 

we  are  already  at  a best  fit).  In  other  words,  we  find  the 

absolute  minimum  of  the  polynomial  in  (2.23)  over  0 ^ 6 £ 1. 

This  search  was  originally  done  by  simply  partitioning  [0,1] 

into  a large  number  of  points  1,000)  and  finding  the  absolute 

minimum  over  that  set  of  points.  This  kind  of  search  proved 

to  be  quite  expensive  computationally,  however,  since  it 

s t 

required  evaluating  an  (N-1)  degree  polynomial  (N 
15-50)  at  more  than  1,000  points,  for  each  iteration  of  the 
fitting  process,  for  each  molecular  species,  in  every  fre- 
quency level.  Therefore,  the  search  has  been  improved 
by  taking  into  account  the  third  derivative  of  the  polynomial. 
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2. 5. 1.1  A Global  Minimization  Procedure  for  Polynomials 

The  problem  is  to  find  the  absolute  minimum  of  P(6) 
for  ® £ 1 • The  method  which  we  shall  give  can,  of 

course,  be  generalized  to  an  arbitrary  interval  [a,b]  , 
and  it  applies  to  any  function  P(6)  which  can  be  decom- 
posed. 


P(0)  = P'^(e)  - P"(6) 


(2.24) 


into  the  difference  of  two  functions  P and  P~  with  mono- 
tone non-decreasing  second  derivatives  on  [a,b]  . An  algo- 
rithm for  performing  this  decomposition  for  a polynomial  is 
given  in  Section  2. 5. 1.2.  Such  a decomposition  is,  of  course, 
not  unique,  since  one  can  add  any  function  with  positive  third 
derivative  to  both  P^  and  P~  . Nevertheless,  because  the 
efficiency  of  the  procedure  depends  on  the  extent  to  which 

d^p' 


a2p+  ^2, 


d02 


exceeds  |p''(0)|  on  [0,1],  not  all  decomposi- 


tions are  equally  useful.  The  one  that  has  been  selected  is 
quite  sophisticated  because  a simpler  one  proved  inadequate 
to  produce  rapid  convergence  to  the  minimum. 


The  algorithm  is  based  on  bounding  properties  of  certain 
quadratic  approximants  to  P'*’(0)  and  P”(0).  Consider  first 
the  quadratic  q'*^(0)  which  matches  P'*’(0)  at  0^^  and  02 
(0  £ 0j^  < 02  £ 1)  and  which  matches  the  derivative  of  P"*^(0) 
at  02  = 


Q'^(0^)  = P'^(0j^)  , 

Q'^(02)  = p‘^(62)  » 
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ar  - 


(e  ) - dP'*'  . . 

de  '®2^ 


de 


By  construction  of  P"*"  , 


d^p"^ 


d0  3 


^ 0 for  0e [0,1] 


The  error  in  the  approximation 


A(0)  - P‘^(0)  - q'^(0) 


therefore,  has  the  properties 


A"  ' (6)  > 0 , 0e[0^,02]  , 


A(0j^)  = A(02)  = A'  (02)  = 0 


Based  on  these  properties,  a proof  is  now  given  that  A > 0 
throughout  the  interval  [62^,82]  . By  Rolle's  theorem,  be- 
cause A(0j^)  = A (02)  = 0 , there  is  a point  ¥e(0j^,02) 
such  that  A'(0)  = 0 . By  Taylor's  theorem,  we  may  expand 
A' (0)  about  0 as  follows: 


A'(0)  = A"(0)(0-0)  + yA' " 1^(0)  ] (0-6)2 


C 


where  ^(Q)  is  a point  between  0 and  0.  Evaluating  this 
2 


at  0 = 0„  , 


0 = A"  (0)  (02-0)  + |a'"U(02)]  (02~0) 
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since  the  second  term  of  this  expression  is  non-negative,  it 
follows  that 


A"  (6)  < 0 

The  case  A' '(6)  =0  is  only  possible  when  P^(6)  is 
quadratic,  for  then  by  Equation  (2.25)  A'(6)  ^ 0 which  is 
irreconcilable  with  A(6j^)  = A(02)  = 0 unless  A (6)  = 0. 

Ignoring  the  trivial  case  when  is  quadratic,  then,  we  have 

A"  (0)  < 0 

so  that  0 is  a local  maximum.  Hence,  there  are  no  local 
minima  in  (0^,02)  and  so  A(0)  must  attain  its  minimum  on 
[01,02]  at  an  endpoint, 

A(0)  ^min[A(0^),  A(02)3  =0  on  I0i»02]  • 

Hence  it  has  been  proven  that  Q^(0)  forms  a lower  bound  for 
P^(0)  on  the  whole  interval  [02f02^' 

q'^(0)  < P'^(0)  , 0e[0^,02l  . 


A similar  proof  will  establish  that  if  the  quadratic  Q~(0) 
is  chosen  such  that 


SSS-R-75-2556 


Q"(02)  = p”(62)  , 


and 


de  de  * 


then  Q (6)  forms  an  upper  bound  for  P”(6)  , 


Q"{6)  ^ P"(0)  , 0el03^,02] 

Clearly,  then,  Q = Q*^  - Q forms  a lower  bound  for  P, 


Q(6)  = q'^(0)  - q"(0)  < p'^(0)  - p”(6)  = P(0) 


on  [0j^,02]  . The  minimization  algorithm  rests  on  this  pro- 


perty. 

Let  us  now  trace  through  a single  iterative  step  of 
the  minimization  algorithm.  Presume  that  the  original  inter- 
val [0,1]  on  which  the  minimum  of  P(0)  is  desired  has  been 
divided  into  subintervals  [0j^,02l,  [03fe4],  t®2N-l'®2N^ 

which  are  the  remaining  candidates  to  contain  the  minimum. 

The  ordering  of  these  intervals  is  such  that  if  Q^(0)  is  thi 
quadratic  approximant  of  the  type  defined  above  for  the  inter- 
val ^®2i-l'®2i^  ' 


qi  = 


min 


0^(0) 


t®2i-l'®2i^ 


then 


— ^2  — ^3  * * * — ^ 


N 
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In  other  words,  the  intervals  are  ordered  according  to  the 
minima  of  their  quadratic  bounds.  Since  the  quadratic  bound 
in  the  first  interval  [92^,62]  dips  lowest,  this  interval  is 
regarded  as  the  most  likely  candidate  to  contain  the  minimum. 

Let  be  the  position  at  which  Q^(0)  attains  its  minimum, 

that  is , 

Qi(yi)  = qi 

Divide  [6^,62]  into  two  new  intervals  and 

establish  quadratic  bounds  separately  for  each  of  the  new 
intervals,  and  insert  the  new  intervals  into  the  candidate 
stack  based  on  their  corresponding  q's.  Drop  the  old  interval 
drop  either  of  the  two  new  intervals  if  its 
quadratic  bound:  (a)  arches  upward  (has  negative  curvature) 
rather  than  dipping  downward,  or  (b)  does  not  have  a local 
minimuiT!!  'Jx-chin  the  interval.  Finally,  drop  any  intervals  k 
for  which 

^ min  P (9 . ) 

l£j<2N  ^ (2.26) 

since  in  such  intervals  P is  bounded  above  an  already-known 
value  of  the  polynomial.  (The  latter  criterion  for  interval- 
dropping  is  particularly  simple  when  the  intervals  are  q-ordered; 
for  if  we  begin  our  search  at  q^^,  and  q^^  is  the  first  q 
such  that  (2.26)  is  satisfied,  then  it  is  also  satisfied  for 

^K+1'  ‘^K+2, 

The  iteration  is  initialized  by  N = 1,  [6j^,02]  = [0,1]. 

It  terminates  when  any  one  of  the  following  convergence  criteria 
is  satisfied; 
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(1)  after  any  interval-dropping  operation,  only 
one  interval  remains  in  the  candidate  stack; 

(2)  more  than  a pre-set  ntunber  of  iterations  has 
been  performed; 

(3)  the  sum  of  the  mimber  of  iterations  performed 
and  the  number  of  remaining  intervals  in  the 
candidate  stack  exceeds  a certain  constant; 

(4)  PCy^^)  < 0 and  q2^/P(y2)  < t some  t > 1 

(we  currently  use  e = 1.01)  . 

The  reasons  for  criteria  (1)  and  (2)  are  fairly  obvious.  The 
third  criterion  was  based  on  the  idea  that  as  the  n\mber  of 
iterations  increases,  fewer  and  fewer  intervals  should  remain 
in  the  candidate  stack  if  the  algorithm  is  functioning  properly. 
The  fourth  criterion  has  been  specialized  to  the  exponential 
fitting  application  because  we  require  the  .polynomial  minimvun 
to  be  negative,  but  the  idea  behind  it,  that  we  are  close  to 
the  minimum  when  the  quadratic  bound  in  the  leading  interval 
closely  approximates  the  polynomial,  is  general.  In  practice, 
surprisingly,  it  is  usually  criterion  (1)  which  terminates 
the  iteration  which  means  that  the  interval-dropping  feature 
is  very  effective. 

Another  version  of  this  algorithm  was  also  developed 
which  required  no  derivative  evaluations  of  P(6)  except  ct 
the  endpoints  6=0  and  0=1,  but  it  proved  to  be  some- 
what less  accurate  and  efficient  than  the  pres«=»nt  version  and 
it  also  necessitated  considerably  more  complex  logic  (especially 
with  regard  to  interval-dropping)  . It  was  based  on  the  fact 
that  the  quadratic  Q(6)  which  matches  any  function  with  non- 
negativt  third  derivative  f(6)  at  three  po  i 
0 £ 01  < 02  < 63  £ 1, 

0(6^)  = f(6.)  (i  = 1,2,3) 
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exhibits  the  "under-over"  property, 

Q(6)  < f(6)  in  [63^,62] 

and 

Q(6)  > f(G)  in  [62,93] 

Hence  by  fitting  p'*'  by  q"*"  at  62  < 63  < 6^  and  p”  by 
Q at  6^<62<63,we  have: 

Q+  _ Q-  ^ p+  _ p-  = p in  162,63]  . 

The  remainder  of  the  algorithm  is  then  similar  to  that  des- 
cribed above . 

Convergence  of  the  algorithm  can  be  proved;  however, 
the  proof  will  not  be  given  here. 


2. 5. 1.2  Splitting  the  Polynomial 

The  crucial  step  in  the  minimization  algorithm  of 
Section  2. 5. 1.1  involves  splitting  P(6)  into  the  difference 
of  two  polynomials: 

P(6)  = P'^(e)  - P"(6) 

with  non-negative  third  derivatives, 

d^P± 

|q~  >0  for  6e[0,l] 

The  algorithm  which  we  use  to  accomplish  this  splitting  will 
be  illustrated  by  actually  decomposing  several  sample  poly- 
nomials. 
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Suppose  that  d(6)  = P* *'(&).  Then  we  want  to  decom- 
pose D{6)  into  a difference: 

D{6)  = D'^{e)  - d"(6) 


such  that 

D-(0)  ^ 0 for  ee[0,l] 

In  order  to  split  D(6)  into  positive  and  negative  (D  ) 

parts,  we  shall  make  use  of  the  elementary  observation  that: 

6^  £ 6^  when  m > n and  6e[0,l]  . (2.27) 

Then  if  we  have,  for  example,  the  pair  of  terms 

36^  - 56' 

we  can  take  the  part  -36'  of  the  negative  term  and  include 
it  with  the  positive  term, 

36'  - 36' 

and  still  have  an  expression  which  is  non-negative  on  0 £ 6 ^ 1. 
For  this  simple  example,  then, 

D'^{6)  = 36’  - 36'  , 

and 

D“{6)  = 26' 

Proceeding  to  a more  complex  example,  consider  the 
quadratic: 

D(6)  =1-6-6^ 
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Divide  D(0)  into  blocks  in  each  of  which  the  coefficients 
are  mono-signed: 

1 - e - 0^ 

block  1 block  2 

The  second  block's  coefficients  are  multiplied  by  a constant 
such  that  their  sum  equals  the  negative  of  the  coefficient 
sum  for  the  first  block: 


This  expression  is  then  non-negative  by  the  property  (2.27), 
and  furthermore  there  are  no  more  coefficient  blocks  to  process. 
Note  that  the  expression  (2.28)  vanishes  at  S = 1,  a property 
whicl\  we  build  into  it;  therefore,  it  may  be  factored, 

(1-  0)  + |e)  . 

Since  the  second  factor  is  strictly  positive,  we  are  finished 
(in  the  next  example  this  will  not  be  so)  . The  splitting  is; 

D^(0)  = 1 - je  - ^0 

and 

D"(0)  = ie  + ie 

An  example  with  an  odd  number  of  mono-sig,,ed  blocks 
introduces  the  further  complications  of;  (a)  keeping  track  of 
'remainders',  and  (b)  performing  a second  blocking  operation 
on  what  is  left  after  the  first  factor  of  (1-0)  is  removed. 
Consider: 

D(0)  = 2 - 79  + 60^ 
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The  mono-signed  blocks  are  merely  single  terms  here.  We  can 
take  the  part  -26  of  the  second  term  and  include  it  with  the 

first  term  to  yield: 

2 - 26  - 56  + 66^ 

Tie  first  two  terms  now  form  a non-negative  expression.  By 
dividing  further, 

2-26-56+56^+6^  , 


the  second  two  terms  can  be  made  into  a non-positive  expres- 
sion. The  extra  0^  which  is  left  over  has  no  terms  to 
match  with  it  and  lienee  is  called  a 'remainder.'  All  remainders 
are  strictly  non-negative  or  non-positive.  They  are  shunted 
off  into  a 'remainder  table'  during  the  blocking-factorization 
process,  and  each  repetition  of  this  process  creates  (in  general) 
a new  remainder  for  the  table.  After  all  possible  blocking- 
factorization  processes  have  been  done,  the  remainders  are 
re-assembled  into  either  d"*^  or  D » depending  on  their  sign- 
As  a shorthand  notation,  we  shall  keep  the  remainders  to  the 
right  of  the  polynomial,  so  for  our  current  example: 

2 - 26  - 56  + 56^  6^ 

Factoring, 

(1  - 6)  (2  - 56)  6^ 


Now  we  block  and  separate  the  second  factor  in  the  same 
fashion: 

(I  - 6)  (2  - 26  - 36)  6^ 
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Shunting  off  the  remainder, 

(1  - 0)  (2  - 20)  QS  -30(1  - 6)  . 

Now  we  are  finished,  since  the  second  factor  can  be  blocked 
and  separated  no  further.  Re-assembling  the  positive  and 
negative  parts, 

= (1  - 0)  (2-20)  + 0^  = 2 - 40  + 30^  , 

and 

d”  = -30(1  - 0)  = -30  + 30^ 

As  a final  example,  consider; 

1 - 20-0^ 

12  3 


The  polynomial  is  separated  into  three  mono-signed  blocks  as 
indicated.  Since  the  largest  negative  value  of  the  second 
block  (at  0=1)  is  -3,  we  can  take  ^ of  the  second  block 
and  adjoin  it  to  the  first  block  and  still  have  a non-negative 
first  block: 


1 2 3 


Since  the  largest  negative  value  of  block  2 is  now  -2,  we  can 
adjoin  j block  3 to  it  and  still  have  a non-positive  second 
block : 


i 
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1 

( 

i 

1 

I 


\ i 

j 

) 

'■  ^ 

’ i 

1 


i 


Block  3 now  has  nothing  to  match  with  it,  and  so  is  placed 
in  the  remainder  table.  (1-0)  is  factored  out  of  blocks 
1 and  2 : 

(1  - 6)  II  + y0  - (|e  + 20^) ] 03 

1 " ^ ^ 

We  now  have  new  blocks  1 and  2,  as  indicated,  and  the  whole 
process  begins  again  (note  that  the  ^0  and  -j0  are  not 
combined  — this  is  in  order  to  avoid  doing  the  actual  factori- 
zation in  the  computational  implementatioi  of  this  method) . 

The  value  of  blocks  1 and  2 at  0 = 1 are  ■j  and  respec- 
tively, so  we  can  adjoin  of  block  2 to  block  1 and  still 
have  a non-negative  first  block: 


(1  - 0)  II  + l0  - 4^(|e  + 202) 


+ 202)] 


02 


Block  2 has  nothing  to  match  with  it,  and  so  is  added  to  the 
remainder  table.  Another  (1  - 0)  is  factored  out  of  block  1, 
to  yield: 


(1  - 0)  (1  + 1^0)  03^  _ |_(i  _ 0)  (4g  ^ 202) 

The  factor  multiplying  (1-0)  is  positive,  so  we  are  finished. 
We  may  therefore  identify: 

D'^  = (1  - 0)'u  + ^0)  + 03 


= 1 - 1.20  - 0.602  + 1.8Q2  , 

and 

^ 29^)  = 0.80  + 0.402  - 1.202  ^ 
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As  noted  above,  the  computational  implementation  of 
this  algorithm  avoids  doing  the  actual  factorizations  in- 
volving (1  - 0)  by  noting  that  what  we  really  want  out  of 
the  factor  multiplying  (1  - 6)^  for  each  n is  the  value 

of  each  of  its  blocks  at  0=1.  This  can  be  simply  re- 
t h 

lated  to  the  n—  derivative  of  blocks  of  the  'reduced' 
polynomial  (the  polynomial  with  remainders  removed)  at  0=1. 
Therefore,  these  derivatives  ere  calculated  instead  of  the 
factorizations,  '^his  saves  a substantial  amount  of  computer 
time . 

2. 5.1.3  Exponential  Fits  of  Representative  Transmission  Functions 

In  order  to  illustrate  the  type  of  exponential  fits 
to  transmission  functions  which  are  generated,  sample  calcula- 
tions are  shown  in  Tables  2.1  to  2.3  for  the  five  frequency 
intervals  180  - 240  cm  ^ (far-infrared  water  vapor  rotation 
band),  720  - 740  cm"^  (CO..  15y  band),  800  - 840  cm"^  (8  - 
12p  "window"),  5440  - 5760  cm  ^ (near-infrared  J^-band  of 
water  vapor) , and  32000  - 33000  cm  ^ (Huggins  band  of  ozone) . 
These  intervals  are  representative  of  the  parts  of  the  spec- 
trum in  which  they  lie,  and  will  serve  to  point  up  peculiari- 
ties of  the  fitting  process.  Table  2.1  contains  the  values 
T^^(nAu),  n = 0,1,2,  ...,  n^,  of  the  transmission  function  to 
which  the  exponential  sum 

E-k^u 
aie 

is  to  be  fitted,  the  corresponding  values  E^^(nAu),  the  per- 
cent difference  between  E.  and  T.  , and  the  coefficients 
a^  and  exponents  k^  resulting  from  the  fitting.  Tables 
2.2  and  2.3  illustrate  the  effects  of  changing  certain  fitting 
paramete.'.o , and  contain  only  coefficients  a^  and  exponents 
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The  range  of  transmission  values,  1.0  - T^^(n  Au) , 
used  in  the  tables  is  not  fixed,  but  depends  on  the  maximum 
estimated  effective  amount  u*^^  of  each  absorber  which 
could  possibly  be  encountered  along  a slant  path  making  an 
angle  of  80“  with  the  zenith.  (Actually,  maximum  H2O, 

CO2,  and  amounts  have  been  taken  to  be  twice  the  current 
values  in  order  that  our  exponential  fits  nught  be  applicable 
to  changed  future  or  past  climatic  conditions.)  Thus,  for 
example,  in  the  800  - 840  cm  ^ case  of  Table  2.1(c),  the 
range  of  water  vapor  transmission  fitted  is  0.775  to  1.0 
because  = 0.775  for  water  vapor  in  this  spectral 

interval.  At  the  other  extreme  are  spectral  intervals  such 
as  180  - 240  cm  (see  Table  2.1(a))  for  which  T,  (u*  ) 

xs  orders  of  magnitude  below  0.001,  the  smallest  transmission 
predictable  with  McClatchey's  scheme  ^^®^which  sets  any  trans- 
mission below  0.001  to  zero.)  Even  should  we  manage  to  ex- 
tend the  transmission  data  below  0.001,  however,  it  would  not 
be  desirable  for  two  reasons  connected  with  the  numerics  of 

First,  because  the  fitting  scheme  takes  equal  steps 
Au  of  absorber  amount  between  the  transmission  data  points, 
fitting  to  transmission  ranges  even  as  large  as  0.001  - 1.0 
results  in  taking  many  data  points  below  0.1  and  relatively 
few  between  1.0  and  0.1,  so  that  the  transmission  function 
is  inadequately  resolved  in  0.1  - l.o  and  the  fit  is  rela- 
tively poor  there.  Secondly,  the  fitting  algorithm  becomes 
unacceptably  poor  if  more  than  about  125  data  points  are 
used,  so  that  increasing  the  resolution  in  0.1  - 1.0  by 
taking  a larger  n^  is  also  unfeasible.  We  have  temporarily 

solved  this  problem  by  putting  a lower  bound  Tr  . (usually 

A AAc  A nt  \ min  •' 

0.005  or  0.01)  on  the  transmission  data  actually  fitted. 

The  variation  of  the  fit  with  Tr  . will  be  discussed  below. 
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Evans*  believes  that  we  can  relax  the  equal  Au 
restriction  by  taking  M steps  of  Au  (starting  ,;t  u = 0)  , 
then  M steps  of  2Au,  then  M steps  of  4Au,  etc.,  where 
M is  the  maximum  number  of  terms  we  expect  the  exponential 
fit  to  contain  (M  = 20  for  the  band  model  we  use  in  the  far 
IR,  M = 8 for  McClatchey's  transmission  data).  This  in- 
volves some  rather  extensive  code  modifications,  however, 
and  has  not  been  done  because  of  more  pressing  problems. 

Once  done,  it  would  eliminate  the  problem  of  adequately  re- 
solving over  its  full  range  of  variation. 

The  number  of  values  n^  of  the  transmission  which 

are  used  for  fitting  is  made  a function  of  the  smallest 

transmission  value  t . = maxfT.  (u*  ) ,Tr  1* 

mm  Av  max''  min-*  * 

n.  =t.  n.  +(l-t.)n 
t min  min  ' min'  max 

where  n^^^^  = 5 for  Tables  2.1-2. 3 and  n^^^^  = 40  for  Table 

2.1  and  80  for  Tables  2. 2-2. 3.  This  leads  to  computational 
savings  when  t^^  is  near  unity  and  assures  that  a full 
”max  Poi^^ts  are  used  when  is  near  zero.  The  linear 

nature  of  this  formula  is  not  optimal,  however,  for  it  leads 
to  unnecessarily  large  values  of  n^  when  • 

Further  study  is  needed. 

Tables  2.1  - 2.3  contain  information  about  both  the 
underlying  line  structure  in  each  spectral  interval  and  about 
the  nature  of  the  fitting  process.  We  begin  by  comparing 
the  various  tables  as  regards  the  fitting  process.  In  the 

spectral  interval  180  - 240  cm  the  transmission 
data  are  supplied  by  a Goody  random  band  model  (which  is  used 
for  V < 340  cm  ^ in  ATRAD) ; for  the  other  spectral  intervals, 
the  tabular  transmission  data  of  McClatchey  are  used.  Because 

* 
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Table  2.2 

Coefficients  and  exponents  for  the  exponential  fits  of  transmission  function 
in  the  five  spectral  intervals,  180  - 240  cm“J-,  720  - 740  cm“l,  800  - 840  cm 
5440  - 5760  cm  and  32000  - 33000  cm  ^ (n  = 80.  Tr  . = 0.04^ 
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of  the  different  origins  of  the  data,  there  are  rather 
striking  differences  between  both  the  accuracy  of  the  fits 
and  the  number  of  terms  in  them:  the  fit  in  Table  2.1(a) 

is  exact  to  the  accuracy  carried  by  our  computer  and  con- 
sists of  18  terms;  the  fits  in  Table  2.1(b)  - (e)  are  con- 
siderably more  approximate  and  consist  at  most  of  6 terms. 

As  a general  principle,  the  larger  tVie  number  of  continuous 
derivatives  possessed  by  the  transmission  function  T^^(u) 
generating  the  data,  the  closer  will  be  the  exponential 
fit.  This  is  obvious  from  an  intuitive  standpoint  because 
a sum  of  exponentials  is  infir itely  differentiable,  and  so 
can  only  match  exactly  with  another  infinitely  differen- 
tiable function  belonging  to  the  function  space  of  exponen- 
tial sums.  The  used  in  the  band  model  is  in  fact 

infinitely  differentiable  on  0 _<  u < “ , while  the  Mc- 
Clatchey  data,  involving  as  it  does  linear  interpolation 
in  a table  of  transmissions,  has  discontinuous  first  deriv- 
atives. (This  suggests  one  way  in  which  the  appropriateness 
of  the  McClatchey  scheme  for  exponential  fitting  might  be 
improved,  which  is  to  use  instead  of  linear  interpolation, 
interT>olation  schemes  of  higher  differentiability  such  as 
cubic  splines.) 

A common  feature  of  every  one  of  the  fits  in  Table 
2.1  is  the  seemingly  random  variation  of  the  sign  of  error. 
This  indicates  that  the  exponential  fit  wanders  above  and 
below  the  transmission  data,  as  it  should  if  it  is  a least 
squares  fit.  Inaccuracies  in  the  fitting  process,  caused 
for  example  by  using  too  many  data  points  (like  n = 150) , 
can  often  be  detected  by  observing  this  sign  pattern.  When 
all  the  errors  are  of  one  sign,  for  example,  computational 
problems  are  definitely  indicated. 
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The  following  approximation  may  be  viewed 


a . e 
1 


-k . u 
1 


as  a Lebesque  quadrature  rule  with  coefficient  a.  repre- 
senting the  fraction  of  the  spectral  interval  Av  over 
which  the  absorption  coefficient  is  roughly  k.  . Thus, 
the  exponential  fits  provide  some  insight  into  the  under- 
lying  distribution  of  line  intensity  in  Av  . In  the 
180  - 240  cm"^  interval  of  Table  2.1(a),  for  example,  one 
could  deduce  a sizable  proportion  ('v.  o.2)  of  strong  absorp- 
tion near  k^  = 17386.1  cm^/g  and  a fairly  uniform  distri- 
bution of  absorption  ranging  all  the  way  from  3580  cmVg 
to  24  cmVg.  For  water  vapor  in  the  720  - 740  cm"^  interval 
(Table  2.1(b)),  we  observe  a preponderant  fraction  (0.64) 
of  very  weak  absorption  (k  'v  0.015  cmVg)  , half  as  much 
('v  0.28)  of  14  times  stronger  absorption  (k  'v  0.208  cmVg)  , 
and  a small  fraction  (0.08)  of  yet  16  times  stronger  ab- 
sorption (k  'v  3.28  cm^/g)  . For  ozone  in  the  32000  - 33000 
cm  interval  (Table  2.1(e)),  the  range  of  k's  is  quite 
small,  only  0.63  (atm-cm)"^  to  0.21  (atm-cm)*^,  indicating 
very  little  line  structure.  Of  course,  the  a's,  and  the 
k's  vary  as  we  change  the  number  of  data  points  n , 
the  lower  limit  , and  the  spacing  Au  , so  the  a's 

and  the  k's  cannot  be  taken  too  literally.  Nevertheless, 
qualitative  features  such  as  the  range  of  the  k.'s  remain 
invariant  as  the  details  of  the  fitting  process  are  varied. 
The  only  exception  to  this  rule  is  when  the  fitting  is 
changed  so  as  to  resolve  parts  of  the  T^^(u)  curve  not 
previously  resolved;  in  this  case,  much  larger  and  much 
smaller  k-values  may  arise  than  were  found  using  the  in- 
adequate resolution. 
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To  show  the  effects  of  changing  the  fitting  param- 
eters, Table  2.2  contains  the  coefficients  and  exponents 
generated  for  all  five  spectral  intervals  when  n„.  ^ is 
increased  to  80  (from  40)  and  Tr 


min 


max 

is  increased  to  0.04 


(from  0.01).  Table  2.3  contains  similar  information  for 
-1  __  . -1 


, CO2  in  720  - 740  cm 


and  ozone  in  32000  - 
Tr 


180  - 240  cm 

-1  ^ 

33000  cm  when  n is  left  at  80  and  Tr  . is  further 

max  min 

increased  to  0.15  (the  other  fits  are  rmchanged  from  Table 

2.2)  . 


The  transmission  data  for  the  180  - 240  cm"^  inter- 
val of  Table  2.2(a)  contain  better  but  still  inadequate 
resolution  in  the  range  0.1  - 1.0  as  compared  to  Tible  2.1 
(a).  The  largest  exponent  in  Table  2.2(a)  a factor  of 
than  in  Table  2.1(a) , in  order  to  account  for 
the  initial  steep  decrease  of  the  transmission  which  had 
been  even  more  poorly  resolved  before.  There  are  many  more 
large  exponents  in  Table  2.2(a)  than  in  Table  2.1(a).  A 
continuation  of  this  trend  is  observable  in  Table  2.3(a), 
in  which  the  largest  exponent  has  increased  a further  factor 
of  three  on  account  of  even  better  resolution  in  0.1  - 1.0. 
Note  that  the  smallest  exponent  also  increases  from  Table 
2.1(a)  to  Table  2.2(a)  to  Table  2.3(a),  which  is  due  to  the 


progressive  increase  in  the  consequent  loss  of 


resolution  in  the  tail  of  the  transmission,  where  the 
small  exponents  predomir.ate . Granted  these  quite  under- 
standable changes  in  the  exponent  range,  however,  the  most 
important  thing  to  notice  about  Tables  2.1(a),  2.2(a), 
and  2.3(a)  is  their  qualitative  similarities;  the  number 
of  exponents  and  their  distribution  between  the  extremes, 
and  the  size  and  general  pattern  of  the  coefficients  exhibit 
regularities  which  are  preserved  as  the  range  and  resolu- 
tion of  the  data  change. 
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Comparing  the  water  vapor  parts  of  Tables  2.1(b)  and 
2.2(b)  for  720  - 740  cm  it  is  app^irent  that  only 

change  due  to  the  better  resolution  in  Table  2.2(d)  is  to 
add  a large  exponent  (84.2)  with  a small  coefficient.  The 
other  coefficients  and  exponents  are  virtually  unchanged. 

The  fit  is  independent  of  Tr^.^  since  the  transmission 
range  is  so  small. 


The  CO2  fit  for  720  - 740  cm"^  (a  region  of  large 
CO2  absorption)  in  Table  2.2(b)  is  very  similar  to  that  .f 
Table  2.1(b)  except  that  the  largest  exponent  has  decreased 
by  a factor  of  three  due  primarily  to  fitting  a single 
extra  transmission  datum  between  0.682  and  1.0  (at  0.776). 
The  comparison  CO2  fit  of  Table  2.3(b),  which  fits  trans- 


mission data  having  even  more  resolution  between  0.7  and 
1.0  (namely  points  0.71,  0.76,  and  0.83)  but  no  points  in 
the  tail  beyond  0.15,  adds  a large  exponent  ten  times  bigger 
than  the  largest  exponent  of  Table  2.2(b),  retains  similar 
exponents  in  the  mid-range,  then  deviates  again  for  small 
exponents.  The  addition  of  the  zero  exponent  with  a non- 
negligible  coefficient  seems  strange,  but  it  is  due  to  the 
fact  that  cutting  off  the  data  at  0.15  allows  the  method 
to  insert  a constant  term  (zero-exponent  term)  into  the  fit 
with  any  coefficient  up  to  0.15  if  it  so  desires,  in  an  at- 
tempt to  fit  the  data.  Without  any  transmission  resolution 
below  0.15,  and  without  some  physical  basis  for  putting  a 
lower  bound  on  the  exponents  (hard  to  come  by  because  of 
lack  of  knowledge  of  line  wings),  there  is  no  sound  reason 


for  rejecting  this  zero  exponent. 
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The  ozone  32000  - 33000  cm  ^ fit  of  Table  2.3(c) 
illustrates  the  danger  of  making  generalizations  about  the 
exponential  fitting  process.  Here  we  have  severely  reduced 
the  resolution  in  the  transmission  tail  by  taking  Tr 
0.15,  and  yet  we  observe  practically  no  impact  on  the^'minimum 
exponent  compared  to  Table  2.2(e).  This  is  quite  different 
from  the  situations  in  Table  2.3(a)  and  (b) , where  knowledge 
of  the  small  exponents  was  reduced  in  the  same  circumstances. 
Ozone  in  this  region  is  of  course  anomalous  in  its  relatively 
small  range  of  absorption  coefficient,  but  it  nevertheless 
serves  to  illustrate  the  point  that  it  is  not  always  neces- 
sary to  consider  very  small  values  of  the  transmission. 

There  ^s  very  little  substantive  difference  between 
the  remainder  of  the  fits  i-  Table  2.2  and  their  analogues 
in  Table  2.1. 

In  conclusion,  we  mention  that  when  lihe  range  of 
transmission  values  which  need  to  be  considered  is  small 
enough  > 0.93  in  the  current  code),  the  full 

exponential  fitting  routine  is  bypassed  in  favor  of  a one- 
term  least-squares  exponential  fit  to  three  data  points, 
which  can  be  done  analytically.  Assuming  that  the  one-term 
exponential  approximation  is 


in  order  to  give  unit  transmission  for  u = 0 , and  defining 


0 -■  e 


-ku*  /2 
max' 


/ 


the  least— squares  residual  becomes 
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u’ 


R = IT.  (-ilfi)  - 6]  + [T,  (u*_)  - e^] 


Av 


L\)'max‘ 


The  condition 


dR 

d0 


- 0 then  yields  a cubic  equation  for  0 , 


u_ 


203  + [1  - 2T,  (u*  )]0  - T.  = 0 

Av  max  ■'  Av  ^ 2 ' ’ 

which  can  be  shown  to  have  only  one  positive  real  root  for 
' s of  interest.  An  example  of  such  a one-term  fit,  for 
CO2  can  be  seen  in  Table  2.1(c). 


2. 5.1.4  Tables  of  Fitting  Parameters 

It  is  desirable  to  make  ATRAD  as  computationally  ef- 
ficient as  possible,  with  a view  to  executing  it  a large 
number  of  times  for:  (a)  heat  budget  studies  involving  a 

large  fraction  of  the  globe,  (b)  diurnal  cycle  studies,  and 
(c)  parameter  studies  involving  the  albedo,  aerosol  density, 
cloud  cover,  sun  angle,  etc.  In  none  of  these  multiple 
executions  would  there  be  any  need  to  chance  the  ATRAD 
spectral  interval  structure  - and  since  the  transmission 
function  fitting  depends  on  the  spectral  intervals,  the 
fitting  parameters  can  be  computed  once  and  for  all  and 
kept  in  tables.  This  is  especially  important  in  view  of 
the  fact  that,  in  the  original  version  of  ATRAD,  the  pacing 
items  as  regards  computer  time  were  the  Mie  calculation 
and  the  fitting  calculation. 

Therefore,  the  fitting  is  now  done  by  a separate  code 
module,  called  EVANS-TABLES . The  cable  currently  used  for  ATRAD 
consists  of  fitting  parameters  for  120  spectral  intervals  spanning 
the  spectral  range  of  60  - 48500  cm”^.  ATRAD  has  the  capability 
of  using  the  full  fitting  tables,  or  any  sub-set  thereof.  Using 
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these  tables,  it  is  possible  to  perform  a complete  ATRAD 
cjlculation  (without  Mie  scattering)  for  40  levels,  and  for 
from  6 to  12  Gaussian  angles,  at  a cost  of  anywhere  between 
300  and  400  CPU  seconds  on  the  UNIVAC  1108.  If  edits  of 
each  spectral  interval  are  suppressed,  this  reduces  to  200  - 
250  CPU  seconds.  Clearly,  the  use  of  fitting  tables  effects 
tremendous  computational  savings  over  the  earlier  version  of 
ATRAD,  which  required  40  - 60  minutes  for  a complete  dear- 
sky  calculation. 


2.5.2  Treatment  of  Overlapping  Bands 


In  many  frequency  groups  Av,  more  than  one  atmospheric 
onstituent  may  contribute  to  the  absorption.  ATRAD  then  fits 
the  transmission  function  T^^^  of  each  individual  constituent 
with  an  exponential  sum. 


(u) 


(n) 


Superscript  1 refers  to  H2O,  superscript  2 to  the  uniformly 
mixed  gases  CO2,  etc.,  and  superscript  3 to  O^.  The  product  of 
these  indiviaual  transmission  functions  is  taken  as  the  total 
transmission. 


T 


Av 


ZEE  ^ 


1 J X. 


i=l  j=l 


=1 


2-58 


SSS-R-75-2556 


Then  M - monochromatic  problems  are  calculated, 

one  for  each  term  of  the  triple  sum,  in  each  of  which  the  opti- 
cal depth  of  zone  Iz,z']  is 


At(z,z')  = (z,z' ) + (z,z') 


+ kP^u^^^  (z,z') 


and  the  line  absorption  coefficient  at  z is 


At  (z,z* ) 


k<i>  3 

1 


3z 


(z,z') 


- ^ 

1 3z 


u^^^  (z,z  ' ) 


v(3) 

‘'I 


3 

3z 


(z,z') 


The  2-derivatives  of  the  u's  are  simply  the  integrands  in  the 
expressions  (2.22)  defining  the  u's. 


2.6  SEPARATION  OF  INTENSITY  INTO  SOLAR  AND  DIFFUSE  PARTS 


Because  the  solar  beam  and  the  specular  reflection  (if 
any)  of  the  solar  beam  are  essentially  6-functions  iu  angle,  it 
is  customary  before  solving  the  radiative  transfer  equation  nu- 
merically to  eliminate  these  beams  from  explicit  consideration 
via  the  following  transformation: 


( T - Tspec 

i = { ^ 

V I T - T^olar 


-1  < \i  < 0 

0 < y < 1 


(2.29) 


i^,  the  diffuse  intensity,  may  be  expected  to  behave  'smoothly' 
as  a function  of  angle,  that  is,  to  be  amenable  to  angular 
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quadrature  schemes  which  assume  i can  be  approximated  by  a 

_ 1 ^ T“ 

low-order  polynomial  in  y . The  solar  intensity  1 and 

the  specularly  reflected  solar  inlensity  Y^pec  obtained 

by  considering  the  incident  solar  beam  to  be  acted  upon  only 
by  extinction  k and  specular  r<f lection  p'  : 

^solar  ^ 


-soec  (2t  (Zo)-T  (z))/p 

jspec  ^ V j I u V 

V 2ti  V , s ' ‘ 


<S  (u+y  0 ) 


where 


T (z)  = f K (z‘  )dz' 

V / V 

•'0 


is  the  optical  depth  measured  from  the  top  of  the  atmosphere. 

If  the  transformation  (2.29)  is  applied  to  the  radia- 
tive transfer  equation  (2.6),  the  boundary  conditions  (2.10) 
and  (2.14),  and  the  flux  equation  (2.7),  the  results  are: 


3i 


e. 


y— ^ + <i  = a'B  + — ^ 

dZ  V V V V 2 


P^(z,y,y' )i^ (z,y* )dy'  +Q^(z,y) 


6 S _ 


-T„(Z)/P, 


4,  , , ,iJ  , -12t  (Z,)-T  (2)]/p, 


i^(0,y)  = 0, 


0 < y < 1 
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(Im|)B„(T  ) + p'  ,(1p1)1^  (z,,  |p|) 


V g " V , s 


/: 


2 I dy'  Ly ' )iy  (^0  ,y'  )y' 

0 ' 


S -T  (Z  ) /U „ 


, -1  £ y < 0 

(2.32) 


F^(z)  = 2tt 


j v\U, 


z,y)dy  + y^jS^  e 


-1  ( z ) /y „ 

V ' ' Q 


■ ^^oS^p;^g(yo)  ^ 


- [2t  (z„ )-t  (z) )/y . 

V 0 V '■''^0 


(2.33) 


These  are  the  equations  which  are  solved  numerically  in  ATRAD 
by  methods  described  in  Section  2.8. 


There  are  two  functions  in  Eqs.  (2.30)  and  (2.32)  which, 
if  they  were  not  smooth  functions  of  angle,  could  destroy  the 
assumed  smoothness  of  i^  as  a function  of  y;  they  are  the 
phase  function  and  the  diffuse  bidirectional  reflectivity 

P"  (j*  The  strongest  peak  in  is  due  to  forward  scattering 
by  aerosol  particles;  therefore,  ATRAD  incorporates  an  approxima- 
tion wherein  this  sharp  peak  is  replaced  by  a 6-function,  leaving  a 
truncated  phase  function  which  is  much  smoother  (see  Section 
2.7.5).  The  same  procedure  can  be  applied  to  a broadened  re- 
flection peak  in  such  as  might  be  caused  by  a slightly 

disturbed  ocean  surface.  That  is,  the  peak  will  be  approx imated 

by  a 6-function  and  lumped  into  the  specular  reflectivity  o'  , 

V ^ s 

leaving  a truncated  and  much  smoother  diffuse  bidirectional  re- 
flectivity p”  ^(trunc)  ^ boundary  condition  (2.32). 
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2.7  SCATTERING  TREATMENT 

2.7.1  Rayleigh  Scattering 

The  Rayleigh  volume  scattering  coefficient, 


V ,R 


8tt  ^ 
3 


(n|-l) 

X-'N^ 

s 


2 


6+3»n\ 

6-7pJ 


_E 

kT 


and  the  Rayleigh  phase  function  for  unpolarized  light 


3(l+p„) 


^R^^s^  4 + 2p 


1-p 

1 + m! 


n 


1+P 


n 


are  taken  from  the  review  article  of  Penndorf The  various 
quantities  in  these  equations  are; 


n^  = index  of  refraction  of  air  at  760  mm  Hg 
and  15°C 


X = wavelength 

Ng  = number  density  of  air  molecules  at  760  mm  Hg 
and  15°C  = 2.54743  x lo ‘ ^ cm" ^ 

= depolarization  factor 

k = Boltzmann's  constant 

p ■=  cosine  of  scattering  angle 
s 


is  calculated  according  to  the  new  empirical  formula  of  Peck 
and  Reeder,  which  supersedes  the  older  Edlen  formula. 

p , which  is  a measure  of  the  anisotropy  of  the  air  molecules, 
is  taken  as  0.0303  following  Gucker,  et.al.  ^ 
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2.7.2  Mie  Scattering,  Single  Sphere 

Mie  scattering  of  a plane  wave  of  wavelength  X 

an  aerosol  particle  of  radius  a is  described  rigorously 

13] 

the  following  infinite  senes: 

00 

a ^ ^ (2n+l)  Re(a^+b  ) ( 

ext  271  / j n n 

n=l 

CO 

(2n+l)  da^l^"  + |b  1=^)  ' 

sea  277  X ^ ' n'  ' n' 

n=l 

^abs  ” '^ext  ^sca 


"l  = 


^2 


00 

n=l 

00 


where 


277a 

= -T- 


complex 
- in 


2 


index  of  refraction 


B = ma 
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'i'n  '^n  - rn  ^^(a)ii^(^) 

\(P)  = P^(U) 

^ U7r^{y)  - a-g2)7T^(p)  . 

The  P^  are  Legendre  polynomials  and  the  are  Ricatti- 

Bessel  functions.  / o , and  a , are  the  extinction, 

scattering,  and  absorption  cross  sections  of  the  particle,  i^^ 
and  i2  describe  the  patterns  of  scattered  intensity  polarized 
perpendicular  and  parallel  to  the  scattering  plane  (determined 
by  incident  and  scattered  beams),  respectively.  For  the  un- 
polarized light  treatment  assumed  in  ATRAD,  the  pattern  of 
scattered  intensity  is  proportional  to  i^^  + 12- 

It  can  be  shown  that  is  related  to  the  integral 

of  i^  + i2  over  all  scattering  angles  ; 

'^sca^®^  = (^)  [i^(a,y^)  + i2(a,yg)]dn^  . (2.38) 


T.he  numerical  evaluation  of  the  infinite  series,  Eqs . 
(2.34)  through  (2.3?)  , while  seemingly  straightforward,  has 
sources  of  error  which  are  related  to  the  accumulation  of 
round-off  error  in  recursively  calculating  ip  , C , ir  , 
snd  . Techniques  to  eliminate  these  errors  have  been 
discussed  by  numerous  authors,  Dave,  ^ Kattawar  and 
Plass,^^^^  and  Denman,  Heller,  and  Pangonis  ^ are  the  pri- 
mary sources  of  the  present  ATRAD  formulation  of  the  series 
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summation,  in  subroutine  MIESCT.  Calculated  results  were  com- 
pared with  the  graphs  of  Kerker^  ^ and  the  tables  of  Deirmend- 
jian^^*^^  and  of  Denman,  et.al.^^^^  Agreement  was  found  in  all 
cases  with  the  published  results. 

The  Mie  scattering  calculation  requires  the  aerosol 
complex  index  of  refraction  m^(z)  as  a function  of  frequency 
and  height.  This  is  the  least  v/ell-known,  for  the  real  atmos- 
phere, of  all  the  required  aerosol  data.  Excellent  m^  data 
are  available  for  water;  we  have  selected  the  measurements  of 
Irvine  and  Pollack , Rusk,  et.al.,^^^^  and  Robertson  and 

Williams  for  inclusion  in  ATRAD  (subroutine  lOR)  . Even 

for  water,  however,  the  uncertainty  in  the  real  part  of  m^ 
is  as  high  as  5 percent  and  in  the  imaginary  part  as  high  as 
20-25  percent.  Somewhat  less  accurate  data  are  available 
for  ice,^^®^  aqueous  solutions  of  NaCl,^^^^  carbonaceous  aero- 
sol, ^ and  quartz,  ^ all  of  which  have  been  included 

in  ATRAD.  The  above  substances  have  each  been  studied  by 
several  investigators,  so  at  least  the  data  have  been  measured 
independently.  No  such  independent  evaluation  exists  for  the 
index  of  refraction  measurements  of  Volz,^^^'^^^  who  studied 
sea  salt,  dust,  Sahara  dust,  and  several  samples  of  water- 
soluble  aerosol  material;  nevertheless,  these  data  have  also 
been  included  in  ATRAD. 

Aerosol  materials,  of  course,  have  highly  variable 
compositions  and  are  subject  to  uncertainties  regarding  their 
chemical  and  physical  states. 

A large  number  of  transmission  spectra  have  been  taken 
for  atmospheric  dust  and  cor  individual  components  of  atmos- 
pheric dust.  ^ These  spectra  indicate  strong  dust 

absorption  in  the  atmospheric  window  region  8-13y,  which  could 
have  a pronounced  effect  on  radiacive  heating.  However, 
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more  data  than  the  transmission  spectrum  of  a substance  are 
needed  to  deduce  the  index  of  refraction.  Even  less  is  known 
about  the  index  of  refraction  of  many  other  inorganic  aerosol 
constituents.  As  for  the  organic  aerosol  fraction,  which  is 
considerable , the  lack  of  data  is  almost  complete. 


2.7.3  Mie  Scattering,  Polydispersion 

In  reality,  the  atmospheric  aercsol  consists  of 
particles  with  a wide  range  of  radii  a.  T1  xs  variation  is 
described  by  a size  distr i bution  function , 


N (z)  n (a, z) 
aer  v / ' 


a . < a < a 

min  — — max 


where  N (z)  is  the  total  number  density  of  aerosol  particles 
at  height  z and  n(a,z)  is  the  normalized  distribution  of 
radii  at  height  z,  defined  such  that 

n(a,z)da  = fraccion  of  particles  at  height  z 
with  radii  in  (a,a+da) 


r 


max 


n (a, z) da  = 1 


min 


The  separation  of  N from  n is  convenient  because  in 

practical  applications  we  often  assume  n is  independent  of 
z . The  volume  scattering  and  absorption  coefficients  for 
an  aerosol  described  by  such  a si;.,e  distribution  are 


3 w(z)  = N (z)  f o (a)  n(a,z)da 

v.M  aer  / sea  ' » / 

•'a  . 
min 


(2.39) 
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Xa 

a , (a)  n(a,z)da  . (2.40) 

abs 

min 

The  pattern  of  scattered  intensity  will  be  proportional  to 

y.  a 

[i  (a,u  ) + i„(a,u  )]  n(a,z)da  . 

IS  z s 

3 • 

min 

Normalizing  this  so  that  its  integral  over  all  scattering  direc- 
tions is  4Tr,  we  obtain  the  Mie  phase  function 


/: 


max 


[ij^(a,Ug)  + i2(a,Pg)]  n(a,z)da 


min 


/ / 

•'47T  •'a 


max 


a . 
min 


[ij_(a,pg)  + i2(a,yg)]  n(a,z)da 


Inverting  the  order  of  integration  in  the  denominator  and 
employing  the  relationships  of  Eqs.  (2.38)  and  (2.39), 


D 

*■  V,M 


a 

max 

min 


[i^(a,yg)  + i2(a,yg)]  n(a,z)da 

(2.41) 


Various  analytic  forms  for  n(a,z)  have  been  proposed. 
Perhaps  the  most  popular  is  the  Junge  power-law  distribution 
(ca”^) ; in  the  American  literature  Deirmend j ian ' s modified 
gamma  distribution  ^ (ca^  e for  haze  and  cloud  is  also 

used  extensively.  Log-normal  distributions 


2-67 


SSS-R-75-2556 


-In^  ( 


a /2lr 


r/r^)/2o^j 


are  used  by  some  European  and  Russian  aerosol  researchers. 

f55  1 

Twomey , et.al.,  employ  Gaussian  distributions,  but  these  are 

not  in  general  use  because  most  available  data  indicate  an 
asymmetrical  distribution  of  aerosol  sizes.  With  the  exception 
of  the  Junge  distribution,  which  is  unrealistically  monotonic 
for  all  a,  the  data  do  not  favor  one  size  distribution  formula 
over  another.  All  of  the  above  are  included  as  options  in 
AT RAD. 

The  dependence  of  ^^er^^^  ^ cloudless 

atmosphere  has  been  investigated  extensively ; ^ ^ 
certain  regularities  have  been  determined.  Eltermann,  ^ for 
example,  has  proposed  a clear  standard  atmosphere  with  a pre- 
scription for  N „(z)  which  has  been  included  in  ATRAD  as 

[53] 

an  option.  In  clouds,  N may  vary  widely,  but  most 

models  ignore  this  and  assume  homogeneous  clouds.  The  varia- 
tion of  N (as  well  as  n)  near  a cloud's  upper  and  lower 
boundaries  may  have  important  consequences  for  atmospheric 
heating,  however,  and  should  not  be  ignored.  Such  variations 
can  be  incorporated  quite  straightforwardly,  within  the  frame- 
work of  the  ATRAD  formalism. 

A further  complication  for  realistic  aerosols  is  that 
the  size  distribution  and  the  index  of  refraction  depend  on 
the  relative  humidity,  Hygroscopic  aerosols  may  begin 

to  take  up  water  at  a relative  humidity  as  small  as  50  percent. 
Hence  it  is  impossible  to  specify  realistic  atmospheric  struc- 
ture data  for  ATRAD  without  taking  this  variable  into  considera- 
tion. 
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2.1. A Numerical  Integration  Over  Aerosol  Size  Distribution 

Techniques  for  performing  the  numerical  integrations 
over  size  distribution  in  Eqs.  (2.39)  - (2.41)  are  far  from 
standard.  Most  investigators  use  trapezoidal  quadrature  on 
the  grounds  that  » *^abs  ' ^1  ' ^2  oscillate  so 

rapidly  and  non-uniformly  with  a that  more  sophisticated 
quadrature  .schemes  are  less  accurate.  This,  of  course,  is  a 
tacit  admission  that  their  integration  meshes  are  not  fine 
enough  to  resolve  the  oscillations.  However,  the  present 
authors  performed  se\eral  integrations  over  size  in  which 
the  oscillations  were  resolved,  and  the  improvement  in  accu- 
racy was  too  small  (never  more  than  a few  percent)  to  warrant 
the  additional  computing  cost  — especially  in  light  of  uncer- 
tainties in  the  Mie  input  data  for  a real  aerosol.  Therefore, 
Romberg  quadrature , ^ which  is  based  on  multiple  trapezoidal 
quadratures,  was  selected  for  use  in  ATRAD. 

[52  ] 

The  comprehensive  study  of  Deirmendjian  proved 

useful  in  ascertaining  suitable  integration  increments  Aa 
for  the  Romberg  quadrature.  However,  Deirmendjian ' s Aa’s 
were  chosen  by  trial-and-error ; he  did  not  determine  a univer- 
sal prescription  for  selecting  Aa's  small  enough  to  ensure 
accuracy  yet  large  enough  to  avoid  wasting  computer  time. 
Studies  by  Dave^®^'^^^  indicate  the  danger  of  selecting  Aa 
too  large  and  shew  that  a Aa  which  is  satisfactory  for  inte- 
grating B and  a,  „ , and  w for  small  angles,  may 

be  unsatisfactory  for  integrating  for  larger  angles. 

However,  Dave,  whose  work  was  restricted  to  water  droplets  in 
the  visible  spectrum,  provides  no  rule  for  automatic  selection 
of  Aa  . 
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Our  Romberg  scheme  for  integration  of  the  fundamental 

Mie  functions  a , ^ , and  i,  + i^  over  aerosol  size 

sea  ext  1 2 

distribution  n(a)  was  abandoned  in  favor  of  a trapezoidal 
scheme  with  a variable  integration  increment  Aa  (a  = particle 

^ 7T3 

radius,  X = wavelength,  a ~ — r — ).  This  decision  was  based  in 
large  part  on  a note  of  Dave's^  in  which  he  pointed  out  the 
computational  economies  of  increasing  Aa  after  the  integra- 
tion covers  a certain  fraction  of  the  particles.  "'n  Dave's 
excimples,  the  integration  increment  was  kept  at  Aa  = 0.1 
until  the  fraction  0.99  of  the  particles  had  been  processed, 
then  increased  to  Aa  = O.T  . This  is,  of  course,  the  sim- 
plest type  of  variable  - Aa  scheme.  A somewhat  more  compli- 
cated variant  is  now  used  in  our  code;  it  will  be  described 
below.  The  former  Romberg  scheme,  while  desirable  for  ob- 
serving the  convergence  of  the  integrals,  was  not  only  com- 
putationally cumbersome,  but  required  fixed  Aa  . Hence, 

While  useful  as  a research  tool,  it  could  not  be  retained  for 
doing  large  numbers  of  more  or  less  routine  Mie  computations. 


Our  present  size  distribution  integration  scheme  errs 
on  the  side  of  caution  in  order  to  do  the  extreme  cases 
(A  -►  0)  correctly,  but  in  so  doing  undoubtedly  takes  too  small 
an  integration  increment  in  other  cases.  Fortunately,  when 

the  Mie  calculation  uses  tables  of  ^sca'  ^ext'  ^1  ^2 

the  integration  increment  is  fixed  by  the  tables;  it  need  not, 
therefore,  be  a source  of  concern,  for  even  a greatly  over- 
conservative Aa  is  of  little  import  tc  the  trivial  amount 
of  computation  necessary  to  produce  the  polydisperse  Mie 
functions  from  these  tables.  When  tables  of  o , etc., 

SCa 

are  not  used,  the  initial  integration  increment  is  taken  as: 


Aa  = 


a 

max 


- a 

o An 


min 


(2.4 


. 5 

.1 
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min 
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where  the  interval  of  integration  is  [a  . , a ] . This 

min'  max 

choice  is  based  on  a conversation  with  Dave  in  which  he  main- 
tained that  0.1  was  the  largest  Act  one  can  safely  use,  but 
that  Aa  must  also  be  small  enough  to  adequately  resolve  the 
interval  [a  • , a ] (i.e.,  to  resol-'e  n(a)).  The  incre- 

mxn  luaX 

r.ent  is  kept  fixed  until  a fraction  f of  the  par“:icles  have 
been  integrated  over,  then  the  increment  is  allowed  to  double 
as  many  as  n^  times  before  the  integration  is  finally  termin- 
ated. The  criterion  lor  increment-doubling  is  that  the  maximum 
relative  change  in  any  quantity  being  integrated,  due  to  the 
previous  integration  step,  be  less  than  6^.  An  acceptable 
set  of  parameters  is 

0.99 


0.001  . 


However,  these  are  definitely  too  conservative  for  most  situ- 
ations. 

The  effects  of  varying  Aa  , f , , and  6^  were 

examined  in  several  spectral  regions  for  the  Arctic  stratus 
cloud  whose  size  distribution  is  shown  in  Figure  3.8.  It  was 
found  that  an  increment  Aa  = 0.1  is  indeed  necessary  to  ac- 
curately integrate  , etc.,  in  the  visible  and  down  to 

X 'v,  0.3y  , that  somewhat  larger  increments  (Aa  = 0.2  at  ly) 
are  permissible  as  the  wavelength  increases  into  the  near  IR, 
but  that  progressively  smaller  increments  (Aa  = 0.05  at  lly) 
must  again  be  used  in  the  IR.  The  ~^^~2o~0  estimate  of 

Equation  (2.42)  takes  care  of  the  IR  quite  well,  and  the  0.1 
estimate  takes  care  of  the  visible.  Nevertheless,  Equation 
(2.42)  is  shown  by  these  studies  to  be  too  conservative  through 
the  near  IR,  a deficiency  which  can  be  remedied  easily. 


and 


f = 


^d  = 
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The  fraction  f must  indeed  be  0.99  for  extreme  cases 
such  as  X 0.3y  - even  values  as  large  as  0.95  lead  to 

error  then.  However,  it  has  been  found  possible  to  decrease 
f somewhat  for  larger  wavelengths.  The  largest  value  of 
6^  which  leads  to  acceptable  accuracy  is  0.01,  and  generally 
the  number  of  doublings  n^  should  be  kept  to  5 or  6. 
Decreasing  6^  to  0.001  or  less  will  usually  suppress  doubling 
for  the  longer  wavelengths  while  allowing  it  to  proceed  for 
the  shorter  ones.  This  is  desirable  since  Mie  computations 
at  longer  wavelengths  are  relatively  inexpensive  anyway. 

The  forward  (diffraction)  peak  in  P for  a par  icle 

V f w 

of  size  parameter  a is  proportional  to 


F(z) 


2 J^(z)j^ 


where  z = asinO  , and  is  the  Bessel  function  of  order 

one.  This  quantity  is  shown  in  Figure  2.3. 

The  actual  forward  peak  will  be  a composite  of  the 

peaks  for  all  ae[a  . ,a  ].  The  sharpest  of  the  component 
^ min  max 

peaks  is  due  to  a = a , and  therefore  in  order  to  resolve 

max 

this  peak,  which  varies  substantially  over  Az  = a A6  = 1 , 
we  pick  an  angular  increment 


A0 

o 


2 a 


max 
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n/t 


Figure  2.3  - Graph  of  F(z) 


Actually,  because  the 


TT 


' s and  T ' s 


entering  the  Mie 


n n 

series  (Eqs.  (2.36)  and  (2.37))  are  pre-calculated  with  a 


piedetermined  angular  mesh  AG 
be  the  closest  multiple  of  A6 


tab 

tab 


to 


A6  is  selected  to 
o 

1/2  a 


max 


also  limited  to  be  less  t’^an  9 (currently  1®) 


m 


A6  is 
The  angu- 


lar mesh  on  which  P,,  „ is  calculated  is  then  ta)cen  as  a 

V , M 


user-supplied  and  even  number  of  steps  of  ^6 
then  4A6 


then  2A6 


o o 

, etc.,  with  the  restriction  that  no  step  is 


larger  than  0 


max 


(currently  2.5®).  The  number  of  steps 


of  each  size  is  even  so  that  a piecewise  Simpson  quadrature 


can  be  used  to  integrate  P 


V,M 


in  case  truncation  (Section 


2.7.5)  or  renormalization  (Section  2.7.7)  is  necessary. 
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2.7.5  Phase  Function  Truncation 

Hansen^  ^ and  Potter ^ have  shown  that  the  forward 

peak  of  the  phase  function  P^  can  be  truncated  without 

significantly  changing  the  results  of  the  transport  calculation. 

Thus,  for  0e[O,O^],  where  6^  is  the  truncation  angle,  the 

real  phase  f\inction  P^(z,0)  is  'replaced'  (to  v/ithin  a 

normalization  factor)  by  a truncated  phase  function  P^  ,,  (z,0) 

in  the  manner  of  Figure  2.4.  The  logarithm  of  P is  taken 

^ V ,M 

to  be  a linear  function  of  6 on  [0,0  ] and  the  values  of 

t ...  4- 

M derivative  are  matched  to  those  cf  P „ at  0 = 0^. 

Neither  Hansen  nor  Potter  propose  a general  method  for 
selecting  0 . Therefore  we  devised  the  following  procedure: 

M 0=0  is  less  than  a prescribed  P (currently 

M truncated.  Otherwise,  truncation  is  attempted 

at  a succession  of  0's  beginning  near  0=0  until  either  a 

successful  truncation  is  achieved  (such  that  P^  ..  (z.O)  < P ) 

t T +-  m 

or  0 exceeds  0 (currently  20°).  It  0^  exceeds  0-^  first, 

^v,M  truncated  at  0*^  without  regard  to  provided 

only  that  P^  j^(z,0)  < P^  z,0).  If  the  latter  procedure  is 

unsuccessful,  P^  is  not  truncated. 

Mathematically,  trunction  amounts  to  approximating  the 
forward  peak,  “ ^v,M'  ^ 6-functinn; 

P = (p  - P^  ) + P^ 

v,M  ' v,M  ^v,M^  ^v,M 

= 4ttF  6 («'  - Tj)  + P^  „ 

v,M 

where  F may  be  determined  by  integrating  both  sides  ov^r  a'l 
scattering  angles 
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F 


)clQ' 


- ]sinG  do 

- P^^j^(z,6)]sin0  d0 


(2.43) 


F is  determined  numerically  (in  subroutine  TRUNCT)  using  a 
Simpson  quadrature. 

Substituting  the  6-function  approximation  into  the 
scattering  terms  of  the  radiative  transfer  equation  (2.6), 
and  eliminating  subscripts  and  unnecessary  arguments  for 
brevity, 


scat,  werms  (^,Q')I(^’)d«'  - 3 I 


(«) 


= I7  y[4TTF  5(Q’-^)  + P^]I(^')dQ’  - 6 I(t2) 
= — 4-/'^  I(ti’)di2'  - B(l-F)Kn) 


(Q , n ' ) I (ti ' ) dQ ' 


P'K"^)  , 


where 


P'  = (1  - F)B  , 
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Therefore,  the  scattering  terms  have  their  oriqinal  form,  but 
P and  3 are  replaced  by  P*  and  3'  defined  above. 

2.7.6  A::imuthal  Integration 

The  aximuthally-averaged  phase  function  . (Eq.  (2*5)), 
not  P^  itself,  is  required  in  the  ATRAD  formulation  of  radia- 
tive transfer.  The  azimuthal  average  of  the  Rayleigh  phase 
function  (Section  2.7.1)  is  computed  analytically  in  ATRAD,  from 

^ ^ f (1  + ^2  Mg)d:J) 

Jq 

= f [l+c^  (UM ' ) ^+2c2Uy ' / (l-y2)  (i-p*  2)cos(}) 

Jq 

+ C2  (1-y^ ) (1-y ' ^)cos^(})] d({) 

= Cj  [1  + C2(yy')"  + I C2(l-y^)  d-y'^)] 

where 

_ _ 3(1  + Pn)  _ _ 1 " Pn 

^1  " 4 + 2p^  ^2  1 + p^ 

The  azimuthal  average  of  the  Mie  phase  function  is 
computed  numerically  using  a Simpson  Quadrature  in  order  to 
reduce  errors  affecting  the  renormalization  procedure  (Sec- 
tion 2.7.7)  . 
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“2.  .1 .1  Renormalizat ion  of:  the  Phasp  Punction 

There  are  two  types  of  renormalization  of  the  Mie  phase 

function  in  ATRAD.  The  first  operates  on  the  Mie  phase  function 

before  azimuthal  averaging,  and  consists  in  multiplying 

all  values  of  P by  the  renormalization  factor 

v,M  ^ 


T2.44) 


where  the  integral  in  the  denominator  is  computed  numerically 
using  Simpson's  Rule. 


The  second  type  of  renormalization  operates  on  the 
azimuthally  averaged  phase  function  P^  and  is  required  by 
the  Grant  and  Hunt  algorithm  (see  Section  2.8)  in  order  to 
preserve  flux  conservation.  In  terms  of  the  Gaussian  angular 
quadrature  mesh  ^ the  numerical  equivalent 

of  the  normalization  condition 

= 2 (2.45) 

is 

m 

Cj  ~ ^ ~ 1/ . . . ^m)  (2-^<^S) 

j = l 

where 
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p..  = p*!  + p^r 

3K  3k  3k 


^Jk  '‘u,M*^'*'j'Pk* 


‘’jk  = ‘‘v.H^^-l'j'-Pk’ 


There  are  several  reasons  why  P , „ will  in  general  fail  to 

V / 

satisfy  Eq.  (2.46).  One  is  that  Gaussian  quadrature  does  not 
accurately  represent  Eq.  (2.45);  that  is,  may  not  closely 

approximate  a polynomial  in  p.  Also,  errors  are  introduced  in 
each  step  of  the  process  by  which  is  constructed 


(a) 

the  numerical  computation 
series  (2 . 34  ) - (2 . 37  ) 

of  the  Mie 

(b) 

the  numerical  integration 

over  the 

aerosol  size  distribution. 

Eq.  (2.41) 

(c) 

the  numerical  integration 

to  obtain 

the  truncation  factor  F, 

Eq.  (2.42) 

(d) 

the  numerical  integration 

over  azimuth. 

Tests  have  bee"^  run  to  determine  ATRAD's  errors  from  each  of 
the  four  sources  (a) -(d);  it  has  been  determined  that  in  no 
case  do  they  exceed  1%.  Hence,  the  necessity  for  renormaliza- 
tion will  be  almost  entirely  due  to  the  inadequacy  of  Eq.  (2.46) 
as  an  approximation  to  Eq.  (2.45). 

Let  us  assume  that  the  corrected  values  satisfy- 

ing Eq.  (2.46),  are  to  be  obtain-^d  from  the  calculated  values 
Pjj^  by  applying  a multiplicative  correction. 


f jk  = '1  ^ "jk> 


calc 

jk 
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Unfortunately,  there  are  ^ ^ ^ unknowns  (the  number  is 

reduced  from  by  the  symmetry  requirement  = tjj^) 

there  are  only  m equations  (Eqs.  (2.46))  to  determine  them. 

This  underspcclf ication  can  be  handled  in  many  ways,  no  one  of 
which  seems  much  preferable  on  physical  grounds  over  any  other. 
Grant*  corrects  only  the  diagonal  elements,  = £j 

so  has  a determinate  set  of  equations  for  the  He  points 

out  that  the  phase  matrices  are  usually  strongly  diagonally 

dominant  and  that  his  procedure  thus  incorporates  the  entire 
correction  where  it  makes  the  smallest  relative  change. 

An  alternative  procedure  has  been  investigated.  Assume 

c -1  = e • + 

1 k- 


Then  Eqs.  (4.37)  become 


m 


j=l 


, calc 


VjK>^j 


= 2 - b,  (k«l, . . . ,m) 


where 


b 


k 


m 


E 

j=i 


calc 

pjk 


and  6.,  is  the  Kronecker  delta.  This  is  a determinate  set  of 
]k 

linear  equations  for  the  e^. 

Both  Grant's  method  and  cur  method  of  renormalization 
have  been  included  as  options  in  ATRAD.  The  effect  of  our 


♦Private  communication. 
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method  is,  in  practice,  similar  to  Grant's;  the  largest  cor- 
rections are  made  to  the  diagonal  elements  of  P • , , • Progressive 
smaller  relative  corrections  result  vv’hen  an  clement  is  farther 
from  the  diagonal.  However,  our  method  offers  the  advantage  of 
spreading  the  correction,  while  Grant's  procedure  concentrates 
it  entirely  on  the  diagonal  element.  In  our  opinion  it  is 
difficult  to  justify  correcting  only  a subset  of  the  p^j^'s. 

Several  test  calculations  have  been  run  in  which  the 
Henyey-Greenstein  phase  function. 


L^Sil 

(1  + g2  - 2gpg) 

was  used  for  ^ , with  g = 0.8  and  g - 0.95  so  that 
large  forward  scattering  peaks  were  being  tested.  Twelve 
Gaussian  angles  (m=b)  were  used.  For  these  examples,  it  was 
found  that  reasonably  large  (up  to  8%  for  g=0 . 8 and  up  to  20% 
for  g=0.95)  renormalization  corrections  were  necessary.  The 
results  point  up  to  the  importance  of  truncation  to  ensure 
that  treated  by  low-order  Gaussian  quadrature. 

This  model,  however,  leads  to  a phase  function  with  a 
substantially  altered  area  under  its  forward  peak,  and  since 
this  forward  peak  was  often  truncated  and  the  scattering  co- 
efficient modified  accordingly,  the  net  impact  of  using  the 
H-G  option  was  often  to  cause  a drastically  different  scatter- 
ing coefficient  to  be  used.  In  order  to  mitigate  this  circum- 
stance, the  parameter  g is  now  chosen  so  that  the  area  A 
under  the  real  phase  function  between  0®  and  D°  (computed  by 
a modified  Simpson's  rule)  equals  the  area  under  the  H-G  ap- 
proximation between  0“  and  D® ; 
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(1  + - 2gp) 


u = cos  D 
o 


This  turns  out  to  reduce  to  a cubic  equation  in  g, 


A(2  - A)g^  - 2(1  - A)  [A  + e(l  - A)]g‘ 


- I2A  + A^  + 4e(l  - A)  ]g  + 2(A  - e)  = 0 


where 


e 5 1 - y 


Studies  of  this  cubic  have  shown  that  for  realistic  ranges  of 
A and  e , it  has  three  real  roots,  one  of  which  is  always 
negative,  one  of  which  always  exceeds  one,  and  one  of  which 
lies  between  zero  and  one.  The  last  root  is  the  desired  one. 

The  current  H-G  scheme  has  lost  one  thing  which  made 
the  earlier  scheme  attractive  - computational  speed.  The  old 
scheme  only  required  the  value  of  the  real  phase  function  at 
0 = 0®;  the  new  one  requires  enough  values  to  resolve  the 
real  phase  function  between  0°  and  D® . As  a result,  the  new 
scheme  is  only  about  a factor  of  5 faster  than  a full  Mie 
computation. 


2.7.8  Total  Phase  Function 

The  total  phase  function  is  obtained  from  Pj 

and  P „ according  to 
V ,M 

7T  _ ®v,R  ^v,M  ^v,M 
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where 


is  the  total  volume  scattering  coefficient, 


2.7.9  Tables 


The  Mie  computation  is  by  far  and  away  the  most  expen- 
sive part  of  ATRAD.  For  example,  for  a single  wavelength 
\ = 0.33y  , for  the  Arctic  stratus  cloud  of  Figure  3.8, 

30  - 60  minutes  of  UNIVAC  1108  time  (depending  on  the  angular 
resolution  in  the  phase  function)  is  required.  The  time 
required  is  progressively  less  at  longer  wavelengths,  but 
the  sum  of  all  these  times  for  the  120  spectral  intervals 
currently  used  by  ATRAD  amounts  to  many,  many  hours.  Clearly, 
ATRAD  could  never  be  a satisfactory  research  tool  if  every 
problem  with  clouds  or  aerosols  took  so  long  to  compute. 

The  only  answer  was  to  split  off  the  entire  Mie  computation 
as  a separate  code  module,  which  then  makes  tables  of  o, 


sea' 


’ext' 


and  to  be  read  by  ATRAD.  Using  these  Mie 


tables,  ATRAD  can  now  be  run  for  a typical  cloud  or  aerosol 
problem  in  (UNIVAC  1108)  times  of  the  order  of  5-10  minutes. 


2.8 


SOLUTION  TO  THE  RADIATI'^N  TRANSFER  EQUATION 


The  method  of  solution  of  the  radiative  transfer  equa- 
tion due  to  Grant  and  Hunt has  been  chosen 
to  solve  the  monochromatic  problem  arising  in  this  section, 
Eqs.  (2.30  through  2.33).  Grant  and  Hunt's  procedure  requires 
no  iteration  in  the  presence  of  scattering,  allows  zones  of 
arbitrary  size,  and  is  computationally  economical.  It  brings 
together  several  lines  of  research  in  radiative  transfer.  On 
the  one  hand,  it  represents,  for  1-D  plane  parallel  problems. 
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the  culmination  of  the  discrete  space  formulation  of  transfer 
theory  associated  with  Preisendorfer and  van  de  Hulst. 

On  the  other  hand,  it  reduces  in  the  limit  of  small  zone  size 
to  the  differential  formulation  of  Rybicki  and  Usher, which 
represented  an  elegant  formulation  of  the  method  of  invariant 
imbedding.  Lastly,  it  improves  upon  the  method (by 

eliminating  the  need  for  zone-centered  intensities)  to  the 
point  where  previous  methods  for  plane  parallel  geometry 

are  now  virtually  obsolete. 

The  elimination  of  scattering  iteration  is  a very  signi- 
ficant feature  of  the  Grant  and  Hunt  method*.  Iteration  methods, 
when  optically  thick  clu  ids  are  present,  would  have  been  pro- 
hibitively expensive,  other  attractive  features  of  the  method 
are: 


(1)  a single  inequality  imposed  on  the 
"primary  layer  thickness"  6Xp  (cf. 

Section  2.8.1)  assures  both  positivity 
of  intensities  and  computational  sta- 
bility; 

(2)  flux  is  conserved  provided  only  that 

the  phase  function  is  correctly 

normalized  (cf.  Section  2.7.7);  and 

(3)  error  estimates  are  available. 

More  generally  the  Grant  and  Hunt  method  possesses  a certain 
elegance,  simplicity,  and  freedom  from  ad  hoc  assumptions  which 
alternative  methods  lack. 

The  basis  of  the  Grant  and  Hunt  method  is  the  inter- 
action principle. 
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a*(T2)  = t(i2,T^)  u'^(T^)  + r(T^,T2)  u'(t2)  + z'^(t^,T2) 
u'(t^)  = r(i2,T^)  + t(T^,T2)  u‘(t2)  + z"(t^,T2) 


(2.47) 


Figure  2.5  illustrates  the  various  quantities  involved 
in  describing  the  interactions  in  the  slab  between  and 

T2  . The  total  optical  depth  coordinate  is  t and  y = cos6 
is  the  cosine  of  the  angle  between  the  ray  and  the  vertical. 
The  u* ' s are  positively  directed  (y  > 0)  intensities,  the 
u 's  negatively  directed  (y  < 0)  intensities  (the  y-dependence 
is  implicit).  The  r's  and  t's  are  integral  operators 
representing  reflection  from  and  transmission  through  the  slab 
[Tj^,T2],  and  the  Z's  represent  effective  internal  sources  in 
[1^,12].  We  now  approximate  the  angular  integrals  in  the 
interaction  principle  by  half-range  Gaussian  or  Radau  quadra- 
tures with  weights  c^  and  angles  0 < y^  < . . . < 1 1 

(y  * 1 for  Radau  quadrature).  (By  using  half-range  quadra- 

T791  . 

ture , an  idea  originally  due  to  Krook,^  •*  in  which  separate 
but  symmetric  quadrature  formulas  are  used  for  0 £ P £ 1 
and  1 ^ y ^ 0 , discontinuities  in  the  intensity  across 
y = 0 (0=90®)  which  may  occur  physically  are  allowed  for.) 

The  r's  and  t's  become  matrix  operators,  and  the  u's 
and  Z's  become  vectors; 


u"(t)  = 


,T2  ,±y;^) 


^(Ti ,T2  ,±ym) 


(2.48) 


where  i^  has  been  written  as  a function  of  t rather  than  z 
Eq.  (2.47)  will  now  be  taken  to  be  in  this  matrix- vector  form. 
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For  the  sake  of  computational  economy,  we  intend  to 
use  only  a relatively  few  (m  £ 7)  angles.  Therefore,  the  de- 
sirability of  solving  for  the  diffuse  intensity  i^  rather 
the  full  intensity  , which  includes  the  solar  flux,  is 

apparent;  if  we  used  , many  more  angles  \i^  would  be 

needed  to  resolve  the  solar  peak,  and  these  could  no  longer 
be  distributed  in  the  optimal  (Gaussian)  way.  We  are,  of 
course,  ignoring  a possible  specular  reflection  peak  in  i 

^ V ’ 

which  has  n£t  been  subtracted  out.  (This  would  be  of  concern 
only  over  fairly  smooth  water  surfaces.)  Tests  have  been  per- 
formed for  a variety  of  situations  which  ensure  that  the  angular 
discretization  is  adequate  to  determine  fluxes  accurately. 

2.8.1  Layer  Composition 

We  now  consider  how  r,  t,  E for  each  zone  are  to  be 
obtained.  Each  zone  is  divided  into  sub-layers  called  pri- 
mary layers,  which  are  optically  thin  enough  that  r,  t,  E for 
each  primary  layer  can  be  obtained  from  the  physical  quanti- 
ties  , etc,  of  the  differential  formulation,  Eq . (2.22) 

Then  the  primary-layer  r's,  t’s,  E's  are  "added  up,"  by  a 
process  we  shall  call  layer  composition,  to  obtain  r’s,  t's, 
and  E's  for  the  entire  zone. 

The  appropriate  primary- layer  quantities,  derived  by 
comparing  the  limits  as  At  = i2  - -►  0 of  Eqs . (2.47) 

with  Eq.  (2.30),  are^^'^^ 

t(T2,T^)  = t(x^,T2)  = I - M’^[l  - P'"'(?)c]at  + 0(At=^) 

r(T2,T^)  = r(x^,T2)  = M'^  P'^'(?)cAx  + O(Ax^)  (2.49) 

E"(x^,X2)  " B-(C)Ax  + O(Ax^) 
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To  0(At),  t(x2»Tj^}  and  t(Tj^,i2)  identical;  however,  they 

may  differ  in  O(At^)  terms.  A similar  remark  applies  to  r. 
The  various  scalars,  matrices,  and  vectors  occurring  in  Bqs . 
(2,49)  are  defined  as 


At  = T . 


■^1  1 ^ 1 


M - [Pj  6..]  , c = [Cj  , 

. [F^(t,p. ,p.)]  = 1F^(T,-Pj,-P.)]  , 

P^'(T)  = [P„(T,p. .-Pj))  = i?^(T,-p. ,P.))  , 


w(t) 


gy(^) 

%(T) 


(2.50) 


b(t,p;  = (1  - «(T)j  b^(T(t))  * ^^^^(T,p) 


b-(t) 


B(t ,±Pj) 


B(t,±pJ 


It  is  worthwhile  to  note  that  the  error  terms  in  Eqs . (2.49) 
will,  in  fact,  be  O(At^)  if  ^ is  chosen  as  the  mid-point 
of  ’ fJ^is  follows  from  Taylor's  theorem. 

In  order  to  guarantee  positive  intensities,  all  ele- 
ments of  the  primary- layer  r-  and  t-matrices  in  Eq . (2.49) 
must  be  positive.  Ignoring  O(Ax^)  terms,  this  reduces  to 
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the  requirement  that  the  diagonal  elements  of  t be  positive, 
which  gives  us  the  following  restriction  on  At  : 


At  < At 


max 


min 

l<i<m 


1 


2 


c.  P CC,y. 
1 


(2 


It  is  also  true  that  this  property  of  positive  elements  in 

r71 1 

r and  t is  preserved  by  layer  composition,*^  ■*  so  that  no 
restrictions  beyond  Eq . (2.51)  are  necessary. 

The  question  of  truncation  err-^r,  that  is,  the  error 
incurred  by  neglecting  O(At^)  terms  in  Eqs,  (2.49),  is  also 
of  concern.  Grant  and  Hunt  have  used  as  a practical  crite- 
rion ^ 


At 


At 


max 


Grant*  has  also  made  some  estimates  of  the  growth  of  trunca- 
tion error  during  the  process  of  layer  composition.  He  finds 
that  if  T,  and  are  the  correct  transmission  and  re- 

flection  matrices  for  a layer  of  thickness  nAr  , and  if 
this  is  built  up  by  composing  primary  layers  of  thickness 
At  , then  the  effect  of  small  truncation  errors  n in 
and  p in  R,  is  to  give  approximate  matrices  t , r 
such  that 

I |T„  ■ I < CnX"(|  IpI  I * I In|  i) 

I |R„  ■ 'nl  I ‘ 2nd  IpI  I * I |n|  I)  * 


* 

Private  communication. 


2-89 


SSS-R-75-2556 


where  C and  X are  defined  such  that 


t < CX‘ 
n 


(0  < X < 1) 


and  11  11  is  the  norm  of  Ref.  [71].  Thus,  for  fixed  At  , 

the  absolute  error  in  R^  and  the  relative  error  in 

grow  about  linearly  with  n . If  we  use  an  appro'^imation  in 

which 

IIpII  + llnll  = 0(At“*b 

and  write  n = H/At  for  fixed  H,  then 


1 1T(H)  - t(H) 11  < KCH  X^  At“ 


1 lR(H)  - r(H) 1 1 < 2KH  Ax“ 


(2.53) 


so  that  the  errors  tend  to  zero  as  At  ^ 0 . Currently,  we 
intend  to  ignore  O(At^)  terms  in  Eqs . (2.49),  so  that 
a = 1 (or  a = 2 if  K is  taken  as  the  mid-point). 


The  layer  composition  formulas  are  derived  as  follows 
Suppose  < T2  < . Write  down  the  interaction  principle 

(2.47)  for  the  layer  [t^,T2)  , and  again  for  the  layer 
[t2,t,1  . Eliminate  u‘^(t'2)  and  u (T2)  from  two  of  these 
equations  using  the  remaining  two.  Put  the  first  two  equa- 
tions in  the  form  of  an  interaction  principle  for  the  layer 
[t^,T3]  . Then  identify  r^^  s 
principle.  The  results  are 
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^13  " ^23  ^32^  ^12  ^23 

^31  ""  ^21  ^12^  ^32  ^21 


where 


^31  “ ^32^  ^21 


t - t Ft" 

^13  ^12^  ^23 


^13  “ ^32^  ^^12  ^23  ^12^ 


23 


^13  “ ^12^  ^^32  ^12  ^23^  ^12 


^ ■ ^12  ^32^ 


-1 


^ " ^32  ^12^  " ^ * ^32^  ^12 


(2.54) 


(2.55) 


In  the  special  case  when  the  layer  [t,  ,'c,.]  is 

• + + + - ^ 
homogeneous . (i.e.,  when  w(t) , P (t) , and  P (t)  do  not  vary 

significantly  across  (t^,Tj])  and  if  the  two  layers  [t, ,12] 

and  [32, Tj]  are  of  equal  thickness  At  , then  formulas  (2,54) 

simplify  somewhat: 


r = r, ^ = r-i  - r„  + t r r t 
1 13  31  o 0000 

t,  E t, , = t,-,  = t r t 

1 13  31  000 

Et  = Z*  = t.  r^(r„  E‘  + Z*)  + Z* 
1 13  00^00  0“^  o 

Z'  = Z'  = t^  r„(r„  Z*  + Z')  + Z' 
1 13  00  00  o o 


(2.56) 
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where 


= (I  - r r ) 

O ^ 0 0"^ 


and 


~ ^12  ■ ^21  ^^23  ^^32 


- ^12  ^21  ^23  ^32 


(2.57) 


y±  = r±  = y± 

o ~ ^12  ^23 


Relations  (2.57)  follow  from  Eqs . (2.49)  by  the  assumptions 
of  homogeneity  and  equal  At’s  . The  single-subscript  nota- 
tion indicates  that  r,  t,  Z depend  on  only  one  argument 
instead  of  two,  that  argument  being  the  corresponding  layer 
thickness.  Defining 


r^^  = r(2^  At) 


etc . 


it  is  possible  to  continue  the  above  process,  called  doubling . 
through  any  homogeneous  region, 


^n+1  ■ '‘^n  ^n  ^n 


^r.4.1  = + t„  r„  r t 

n+1  n n n n n 


Cl  - ^n  ^n(^.  ^ ^ 


(2.5S) 


= i.,  r (r  I + l')  I 

n*fl  n n'-n  n n'  ^ 


'n+1  n n n 
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Doubling  is  much  faster  than  simple  layer  addition  as  embodied 

N 

in  Eqs . (2.54).  To  add  2 primary  layers  comprising  a homo- 
geneous layer  requires  2^-1  layer  additions  versus  only  N 
doublings . 

In  the  future  the  code  will  allow  for  arbitrarily 
larger  zones,  which  may  or  may  not  be  homogeneous.  If  a zone 
is  homogeneous,  doubling  will  be  used;  if  it  is  not,  the  code 
will  partition  it  into  sub-zones  which  are  homogeneous,  double 
within  each  of  the  sub-zones,  and  add  the  sub-zones  using 
Eqs.  (2.54),  For  now,  however,  the  code  will  be  constructed 
on  the  assumption  that  each  zone  is  homogeneous.  This  is 
fairly  restrictive,  since  it  implies  that  neither  the  scatter- 
ing nor  the  absorption  coefficient  vary  substantially  across 
a zone.  No  more  than  40  zones  should  be  necessary,  however, 
and  this  does  not  seem  overly  burdensome  (Grant  and  Hunt  run 
with  about  40  zones*) . 

2.8.2  Source  Doubling 

The  nroblem  of  an  inhomogeneous  source  Z~  in  a zone 
which  is  otherwise  homogeneous  requires  special  attention. 

An  inhomoj  jneous  source  region  will  occur  frequently  on  ac- 
count  of  the  factor  e ' ^sun  in  B(t,vi)  and  on  account  of 
the  rapid  variation  of  B^(T)  with  T in  many  spectral  re- 
gions. It  turns  out,  however,  that  the  most  frequently  occur- 
ing  types  of  source  inhomogeneity  still  admit  a modified  form 
of  doubling.  Inhomogeneous  sources  can  be  doubled  right  along 
with  r and  t with  the  resulting  saving  of  computing  time. 

Let  us  begin  by  developing  a general  formula  for 
source  composition  when  the  source  variation  across  an  other- 
wise homogeneous  zone  is  described  by  f(r)  . Referring  to 

* 
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Figure  2.6  and  l;qs.  (2.49),  let  the  source  vector  for  any  one 
of  the  primary  layers  be 


:i  = If  f(q) 


i 1 , . . . , 2 


N 


(2.59) 


+ 

where  Zp  are  independent  of  t , and  is  the  mid-point 

of  primary  layer  i , 


5i  = ,0  ^ iAx  - §1 


(2.60) 


Adding  the  sources  in  layers  1 and  2 according  to  Eq._(2.54), 
it  is  clear  that  the  source  + 2 combined  layer 

will  be  a linear  combination  of  f(C^)  f(^2^* 


i 


........  A...,.-.,-  ...1 


Figure  2.6  — Homogeneous  zone  composed  of  2^ 
primary  layers. 
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f(C^)  + 


f(C2) 


(The  expressions  for  the  Vs  are  not  of  interest  here.) 
Similarly,  the  source  for  the  combined  layer  3+4  involves  only 
a change  in  the  ^'s: 


fCCj)  * £(54) 

Then,  upon  combining  layers  1+2  and  3+4,  a linear  combination 
of  f(Cj^)  through  fCC^)  results: 

4 

l~  = ^ ''  ff£  ) 

^l  + 2 + 3+4  Z-f  SL 

Jl=l 

Suppose  that  any  2^  adjacent  primary  layers  have  been  super- 
posed in  the  above  manner.  Using  the  notation 

^(i,n)  ~ ^i+(i+l)+  +(i+2^-l) 

where  i denotes  the  first  primary  layer  in  the  grouping, 
it  is  possible  to  prove  by  induction  that 

l=l 


We  now  'pecialize  our  considerations  to  the  solar 
source  (f(T)  = 0)  and  the  Planck  source 

{f(t)  = B^(T(t))}  . By  the  linearity  of  the  source  composi- 
tion formulas  (2.54)  it  is  possible  to  compose  the  solar 
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source  and  Planck  source  separately  and  then  simply  add  tlie 
results  at  the  end  to  obtain  the  full  source.. 

For  the  case  of  the  solar  source,  the  quantities 
of  Eq.  (2.59)  are^  from  Eqs . (2.49)  and  (2.50), 


+ S oj(t) 

^P,sol  = -V- 


P^(T,±y^,Uo) 


[_P  ,T,±y  ,Po)  J 


(2.62) 


Since  to  and  are  assumed  not  to  vary  across  the  zone , 

the  argument  t refers  to  the  mid-point  of  the  entire  zone 
in  Figure  2.6.  Now  assume  that,  beginning  with  the  EZ  , 

r , SOi 

and  Eq.  (2.59),  the  first  2 primary  layers  have  been  com- 
posed, so  that  the  E^  are  known.  The  adjacent  2^  pri- 
mary layers  will  then,  by  Eq.  (2.61),  have  source  vectors 


.n 


(l*2".n)  '5^1 


±(n) 


2“  . 


-2"At/u  o . 

= e 


^a,n) 


(2.63) 


n 


^U,n) 
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where 


-2  At/mo 


Remembering  the  definitions  r = r(2  Ax)  , etc.,  the  sources 
^(1  n+l)  combined  layer,  containing  2“  primary 

layers , are 

^(l,n+l)  " ^n  ^n(^n  ^(1+2^, n)  * ^(l,n))  ^ ^(1+2^, n) 

" ^n  ^ni’^n  ^n  ^(l,n)  ^(l,n))  ^n  ^(l,n) 

(2.6 

^(l,n+l)  “ ^n  ^n(^n  ^(l,n)  * ^(l+2^,n))  ^ ^(l,n) 

“ ^n  ^ni’^n  ^(l,n)  ^ \ ^(l,n))  ^ ^(l,n) 


These  are  adaptations  of  the  source  composition  formulas  in 
Eqs . (2.54).  From  the  definition  of  h^  , it  is  obvious  that 


h = h^ 
n+l  n 


Together,  Eqs.  (2.65)  and  (2.66)  form  a doubling  scheme  for 
the  solar  (or  any  exponential)  source.  The  scheme  is  initia- 
lized by 


+ + 

^(1,0)  " ^P,sol  ® 


-(To  T 


ho  =■  e 


-At/p, 
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and  iterated  from  n=0  to  n=N-l. 

In  order  to  derive  doubling  formulas  for  the  Planck 
source,  it  is  necessary  to  have  a more  tractable  functional 
form  for  f(i)  than  B^(T(t))  . It  can  be  shown  that  it 
is  reasonable  to  assume  is  piecewise  linear  in  t , so 

that 


= B.,(T(t))  = B + 


(2.68) 


across  the  zone  of  Figure  2. 6 » where 


®o  " Bv(T(t^)) 


B'  = 


"l  ■ "c 


Expressing  this  in  terms  of  (see  Eq.  (2.60)), 


f(?i)  = * B'(t„  - i/iT  - |i  - T„) 


B^  + B'  iAx 


(2.69) 


where 


The  sources 
come 


v± 

(l,n) 


B'  = B - B' 
0 o 

for  the  first 
2n 


At 

2 


2" 


A=1 
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while  for  the  adjacent 


,N 


primary  layers 


2" 

£=1 


' ^(l,n) 


+ 2”  B'  At  Y" 
n 


i 


(2.70) 


where 


2*^ 


jt=l 


The  Y*  refer  to  a constant  source  f(T)  = 1 , and  so  may  be 
calculated  by  doubling  as  in  Eqs . (2.5-8).  A further  simpli- 
fication is  that  symmetry  prevails, 


Y 


n 


This  is  intuitively  obvious,  since  a constant  source  would  be 
expected  to  emit  the  same  at  -y  as  at  +y  . It  can  be 
proven  inductively  by  first  noting  that  the  Zp  of  Eq . (2.59), 
which  for  the  Planck  source  are. 


^P.plk  ' [1  ■ “(T)]*-'  M 
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are  equal, 

E =T.^  - y~ 

P,plk  " P,plk  “ ^P,plk 

Because  £(t)  = 1 , 


Y“  = Z 
o P,plk 


(2.72) 


• (2.59);  so  that  symmetry  is  proved  for  n = 0 . 
symmetry  has  been  proved  for  a general  n , then  from  the 
doubling  formulas  (2.58) 


"'n.l  = ^n  ^n^^-n  "n  ^ ^n^  ^ "n 


‘n+1 


= t. 


n 


r (r  Y 
n'“  n n 


Y ) 
n-' 


+ Y. 


n 


n+1 


which  completes  the  induction  proof,  and  gives  us  the  simpli- 
fied doubling  formula 


n + 1 


MI  * r^(i 


n-^  ■'  n 


(2.73) 


for  Y . 
n 

Tlie  two  layers,  each  containing  2^  primary  layers, 
are  now  composed  just  as  for  the  exponential  source  (see 
£q_-  (2.65)): 
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■(l,n  + l)  = ^ B'  At  Y^)  + ^(l,n)] 


^ ^(l.n3  ^ 2^  B-  At 


(1 


,n^l3  = "n  ^U.n)  " ^'(l,n)  " B'  At  yJ  . 


(2.74) 

) 


Together,  Eqs . (2.73J-  and  (2.74)  form  a doubling  scheme  for 
the  Planck  (,a:.d  any  linear-in- t)  source.  The  scheme  is 
initialized  by  Eq.  (2.72)\  and  by 


^(1,0)  = * «’  (2-75) 

and  iterated  from  n*0  to  n=N-l, 

More  accurate  piecewise  polynomial  approximations  to 
possible,  at  the  expense  of  more  complicated 
source  doubling  formulas.  In  particular,  economical  piece- 
wise  quadratic  and  piecewise  cubic  schemes  have  been  derived 
and  will  be  programmed  at  a future  date  to  ascertain  if  the 
increased  accuracy  of  the  heating-rate  calculations  so  ob- 
tained warrants  the  additional  computational  c*xpense. 

2.8.3  Forward  and  Backward  Pass e s 

By  the  procedures  of  Sections  2.8.1  and  2.8.2,  we 
obtain  reflection  and  transmission  matrices  and  source  vectors 
for  every  zone  of  the  atmosphere.  If  the  zone  structure  is 


and  the  notation  is 
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^.n±l  = 

^n,n±l  “ 

^n,n+l 


(2.76) 


then  according  to  the  Grant  and  Hunt  algorithm  (whose  deriva- 
tion is  lengthy  and  will  not  be  repeated  here)  the  solution 
for  the  intensity  vectors  splits  into  two  parts. 

The  first  part  is  the  forward  pass . In  the  forward 

pass,  the  m^m  matrices  E , G , and  H and  the  m-vectors 
(1)  (2)  fSI  n ’ n ’ n 

, and  •'  are  calculated  recursively  as  follows: 

= 0 = 0 (2.77a) 


t 

n+1 

■ ^n+l,n 

(2.77b) 

<\< 

n+1 

^n+1  ^n,n+l 

(2.77c) 

n+1 

E H 
n n+1 

(2.77d) 

n+1 

^n,n+l  * ^n+l,n  ^n+1 

(2.77e) 

CD  = 

n+1 

r (r  yC^)  + r"  1 

n+1  ^ n+l,n  ''n  ^njn  + l-' 

(2.77f) 
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V(2) 

n+1 


= + E (2.77g) 

n n n + 1 

' Vl.n  C2.77h) 


(n  = 1,2,  ...,  N) 


The  initialization  of  at  zero  reflects  the  fact  that 

there  is  no  incident  diffuse  radiation  at  the  top  of  the  at- 
mosphere, T = 0 (the  incident  solar  radiation  is  already  ac- 
counted for  by  the  transformation  from  to  i^)  . Economy 

of  storage  is  achieved  by  performing  the  sub-layer  superposi- 
tions, which  yield  the  r's,  t's,  and  Z's,  simultaneously 
with  the  forward  pass;  in  this  way  the  r's,  t's,  and  Z's 
need  not  be  saved  as  a function  of  zone. 

In  order  to  identify  certain  quantities  in  the  back- 
ward pass  formulas,  it  is  necessary  to  look  at  the  surface 
boundary  condition  as  stated  in  Section  2.4. 


^v^^N+1 


>y) 


+ 2 /*  y'  p^(y'  , |y|) 

•'0 

(2.78) 

‘ \^^N+1'^'^  (-1  1 < 0) 


where 


fu(yJ  = e^,(y)  ,(TJ  + -^  y 


V'  g- 


Applying  our  half-range  Gaussian  quadrature  to  the  integral 
in  Eq.  (2.78)  , 
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m 


j = i 


‘‘j  "vf^Nn.Xj’ 


(i  “ 1 , • • • > ni) 

This  can  be  written  in  matrix-vector  form,  using  the  notation 
of  Eq.  (2.76)  as 


where 


“n+1  ^N+1 


w = 


'fbC'p 


LW 


(2.79) 


The  backward  pass  may  now  be  written 


^N+1  " * ^M+1  ^ * ^m4.i  ) 


and 


^N+1  *G 


N+1 


N+1' 


Vl  = ^G 


= G T u'  . + 
n n+1  n+1  n+1 


(2.80) 


u =H  ,u 
n n+1  n+1  ''n+1 


(2.81) 


n = N,N-1 , . . . , 1 
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From  the  components  ’^n  ’ fluxes  are 

obtained : 


»/: 


* “o  S e 


‘ In  ■ >'iv("n-'''Udu 


+ y S e 

^ O V 


(2.82) 


m 

2tt  / , c.  y.[i  (t  ,y.)  - i (t  ,-y  )l 
1 i'-  n’  n’ 

i = l 


This  completes  the  algorithm. 


2.8.4  Special  Case  of  No  Scattering  in  the  Grant-Hunt  Algorithm 

In  the  case  of  negligible  scattering  it  is  not  necessary 
to  build  up  the  solution  from  optically  thin  layers;  simple  ana- 
lytic expressions  for  the  reflection  and  transmission  matrices 
and  Planck  source  vectors  can  be  derived  for  zones  of  arbitrary 
optical  depth.  These  analytic  expressions  are  obtained  from  the 
Grant-Hunt  algorithm  by  taking  the  limit  as  the  primary-layer 

thickness  At  goes  to  zero.  For  the  reflection  and  transmission 
P 

matrices  of  a zone  of  optical  thickness  At  , we  obtain 
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^ij  = 


t,,  = 

in  i-D 


For  a Planck  source  which  varies  linearly  in  x, 


where 


g,  ^ E^ITq-^At)  - B^(t„) 


At 


th 

we  obtain  for  the  j component  of  each  source  vector, 


E 


-Ax/y . 

- (1  - e 3) 


+ 


y .B’ 
D V 


-Ax/y . 

(1  - e 


'plkj  “ (1  - e 


-Ax/y . 


')  + yjB; 


Ax 

1 - (1  + e 

D 


These  special  case  results  are  incorporated  in  ATRAD. 


2.9  FLOW  DIAGRAM  OF  ATRAD 


Figure  2.7  shows  a brief  flow  diagram  of  the  ATRAD  Code. 
EVANS-TABLES  makes  tables  of  the  exponential  fits  of  transmission 
functions.  1112-TABLES  makes  tables  of  the  Mie  scattering  functions 
for  single  homogeneous  spheres  of  varying  sizes,  but  fixed  index  of 
refraction.  MIE-TABLES  makes  tables  of  Mie  scattering  functions  for 
spherical  polydispersions.  The  dotted  line  correcting  1112-TABLES 
and  MIE-TABLES  is  meant  to  indicate  that  MIE-TABLES  has  the  option 
of  either  generating  its  own  values  of  ^ext'  ^1^^2 

reading  them  from  tables  generated  by  1112-TABLES. 
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Figure  2.7  - AT RAD  code  organization,  showing  the  role  of  the 
three  auxiliary  table-making  codes. 
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3 . ATRAD  APPLICATIONS 


This  chapter  presents  the  results  of  comparisons  be- 
tween the  ATRAD  code  and  the  radiation  subroutine  of  the 
Mintz-Arakawa  global  circulation  code.  The  Mintz-Arakawa 
treatment  is  discussed  in  Section  3.1.  The  test  problem  re- 
sults are  presented  in  Sections  3.2  - 3.5. 

3.1  THE  MINTZ-ARAKAWA  RADIATION  SUBROUTINE 

The  calculation  of  radiation  heating  and  cooling  in 
the  Mintz-Arakawa  (M/A)  global  circulation  model  (GCM)  is 
necessarily  a considerably  simplified  one  in  order  to  keep 
the  computational  work  within  practicality.  In  order  to 
determine  the  size  of  the  errors  which  result  from  this  sim- 
plified radiative  treatment,  a series  of  comparisons  of 
fluxes  and  heating  rates  calculated  with  the  first-principles 
radiative  transfer  code,  ATRAD,  and  those  from  the  Mintz- 
Arakawa  two- level  GCM  are  being  made.  The  results  of  the 
first  two  comparison  cases  corresponding  to  cloud- free  con- 
ditions over  land  are  reported  in  Section  3.2.  In  order  to 
make  these  comparisons  we  performed  calculations  of  radiative 
quantities  with  the  M/A  radiative  subroutine  corresponding 
to  a particular  epoch  in  a standard  global  circulation  calcu- 
lation. (We  are  indebted  to  A.  B.  Kahle  of  RAND  for  the 
M/A  computer  code  and  for  input  data  for  it.)  In  the 
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following  section  we  outline  the  features  of  the  M/A  radiative 
model  and  describe  how  the  input  data  for  it  are  transcribed 
to  the  ATRAD  code. 

The  M/A  treatment  as  employed  in  the  climate  dynamics 
program  of  the  RAND  Corporation  is  described  by  Gates, 

[ C ’ J 

et.al.  , who  give  a definitive  description  of  the  formula- 
tion and  code.  V7e  outline  the  principal  features  of  the  code 
here;  howe^'er,  we  are  particularly  interested  in  specifying 
how  the  M/A  atmosphere  is  constituted  in  order  that  we  may 
provide  a proper  initialization  for  ATRAD. 

The  M/A  radiative  subroutine  determines  the  radiative 
heating  rate  in  two  layers  of  the  atmosphere  and  the  radiative 
flux  into  the  surface.  Each  of  the  layers  contains  one-half 
of  the  mass  of  t atmosphere  lying  below  the  200  mb  level. 

The  part  of  the  atmosphere  above  200  mb  is  not  explicitly 
contained  in  the  M/A  calculation  but  is  included  in  ATRAD. 

The  relevant  radiative  quantities  to  determine  atmospheric 
heating  are  completely  specified  by  giving  the  net  radiative 
fluxes  at  the  boundaries  of  these  levels;  in  terms  of  the 
variable  a denoting  the  fraction  of  the  atmosphere  below 
200  mb  these  are  F(a  = 0)  at  the  200  mb  level,  F(a  = 0.5) 
in  the  middle  of  the  atmosphere,  and  F(a  = 1.0)  at  the  lower 
boundary  of  the  atmosphere.  These  three  altitudes  are  in- 
cluded among  the  zone  boundaries  of  ATRAD  for  which  radia- 
tive fluxes  are  calculated. 

The  radiative  model  of  the  M/A  code  distinguishes 
three  spectral  regions  in  which  different  radiative  calcula- 
tions are  performed.  Taking  advantage  of  the  fact  that  there 
is  little  overlap  between  the  solar  and  thermal  radiative 
sources,  the  major  frequency  division  is  between  long-wave 
and  short-wave  radiation.  The  frequency  dividing  these  re- 
gions is  not  specified  in  the  M/A  code  but  is  chosen  to  be  at 
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4y  for  the  purposes  of  editing  data  from  ATRAD.  In  the  long- 
wave region  the  M/A  code  takes  account  only  of  thermal  emis- 
sion and  absorption  by  the  atmosphere  and  the  ground.  In  the 
short-wave  case  where  thermal  emission  is  neglected  the  spec- 
trum is  divided  at  0.9y  into  two  parts;  in  the  near  IR  region 
the  solar  radiation  is  assumed  to  be  absorbed  out  unscattered, 
while  in  the  remaining  spectral  interval  it  is  assumed  to  be 
unabsorbed  but  Rayleigh  scattered.  Using  these  spectral  in- 
tervals, the  results  of  the  M/A  calculation  can  be  presented 
as  a matrix  of  nine  net  fluxes,  corresponding  to  the  three 
spectral  intervals  and  the  three  altitudes. 

The  most  important  atmospheric  absorber  for  typical 
atmospheric  conditions  is  water  vapor.  The  M/A  calculation 
takes  account  principally  of  the  effect  of  water  vapor,  in- 
corporating empirical  corrections  for  the  effect  of  CO2. 

Based  on  a single  value  of  the  water  vapor  mixing  ratio 

at  a = 3/4,  and  an  assumed  saturated  vapor  pressure  above 
120  mb  in  the  stratosphere,  the  code  interpolates  a profile 
of  water  vapor  density.  The  vertical  profile,  Q,  as  a func- 
tion of  pressure,  p,  in  mb  is  given  by 


p > 120  mb  , 


2.06255  X iQ-** 
P 


p < 120  mb 


(3. 


p^  is  the  pressure  at  (j  = 3/4  and  the  constant  k is  deter- 
mined by  ma-ching  the  values  of  Q at  120  mb: 
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ln(Q2/1.7188  X 10"^) 
ln(p2/120  mb) 


Using  this  profile  the  effective  water  vapor  amount, 
including  a factor  for  the  pressure  broadening  of  the  absorp- 
tion lines,  IS  calculated  between  altitude  levels.  This 
quantity  is  subsequently  used  to  form  the  absorptivity  of  each 
layer  for  the  near  IR  short-wave  radiation  and  the  transmis- 
sivity of  each  layer  for  the  IR  radiation. 

Equation  (3.1)  is  also  used  to  form  the  water  vapor 
density  input  data  for  the  ATRAD  code.  In  addition,  it  is 
necessary  to  interpolate  and  extrapolate  the  temperature. 

The  M/A  code  makes  use  of  the  calculated  values  of  temperature 
in  the  upper  layer  and  in  the  lower  layer.  From  these 
values  the  temperatures  at  the  200  mb  level  and  at  the  boundary 
between  the  upper  and  lower  layer  are  found  by  invoking  the 
following  functional  dependence  of  temperature  on  pressure: 


T 


(3.2) 


where  < 


R/C  and 


e = T/p' 


is  proportional  to  potential  temperature.  Equation  (3.2)  in- 
sures that  the  prescribed  values  of  temperature  at  p 

(a  - 1/4)  and  p^  (a  = 3/4)  are  recovered  and  in  addition, 
the  temperature  becomes  zero  at  p = 0 . In  addition,  the  air 
temperature  near  the  ground  at  the  top  of  the  boundary  layer, 
T^,  is  also  calculated  in  the  M/A  code  from  a consideration  of 

3-4 


SSS-R-75-2556 


the  fluxes  of  static  energy  between  the  ground  and  level  3 
(a  = 3/4) . Subsequent  adjustments  of  the  original  input  tem- 
peratures may  take  place  from  convection.  In  order  to  specify 
the  initial  temperature  profile  for  the  ATRAD  code,  we  use 
Eq.  (3.2)  from  0=0  to  a = 3/4  . To  determine  temperatures 
from  0=0. 75  to  o=1.0,  we  use  the  same  expression  ad- 
justed to  match  and  T^: 


T 


e_p . 

3^*1 


84P3 


- 


< 

P', 


P3  1 P i P4  • 


Finally,  in  the  region  above  200  mb,  where  the  M/A  code  gives 

no  guidance,  we  provide  a rough  simulation  of  the  stratospheric 

temperature  rise,  using  an  expression  which  is  continuous 

with  at  p = 200  mb  . 

o 


T = 


T 


o 


tPp 

(p  + kj^) 


0 1 P 1 Pq 


(3.4) 


Having  specified  the  relation  between  temperature  and 
pressure  throughout  the  atmosphere,  we  can  proceed  to  the  deter 
mination  of  the  altitude  by  using  the  hydrostatic  relationship. 


where  is  the  altitude  of  the  surface  and  we  have  neglected 

the  effect  of  water  vapor  in  the  atmosphere.  In  view  of  the 
several  analytic  forms  given  above  for  the  temperature-pressure 
relationship,  we  have  approximated  it  by  linear  expressions, 

T = A + Bp  , within  each  of  the  ATRAD  zones  which  are  of  the 
order  of  25  mb  thick.  The  altitude  of  zone  I becomes 


3-5* 


SSS-R-75-2556 


Z(I)  = Z(I-l)  + |[a  In  + b(pU-1)  - P(I))]  , (3.5 

where 

T(I-l)  P(I)  - T(I)  P(I-l) 

^ Pv'I)  - P(I-l) 

and 

_ T(I)  - T(I-l) 

® P(I)  - P(I-l)  • 

We  have  chosen  Z(l)  = 0,  thereby  measuring  the  altitude  rela- 
tive to  the  surface. 

3.2  CALCULATIONS 

Two  different  locations  on  the  Earth's  surface  were 
chosen  for  comparison  calculations.  The  criteria  for  choosing 
these  locations  were  first,  that  they  be  the  centers  of  grid 
squares  in  the  M/A  code,  second,  that  the  M/A  code  predict 
no  clouds  for  the  chosen  points,  and  third,  that  the  chosen 
points  have  substantially  different  water  vapor  profiles. 

Two  such  locations  were  found  in  the  M/A  run  that  was  furnished 
to  us.  The  first  (Problem  1)  was  at  14°N,  20°E,  in  Chad  (in 
the  Sahara  desert)  . The  second  (P'^oblem  2)  was  at  i8°S,  65®W, 
in  Bolivia.  Temperature  profiles  and  water  vapor  profiles 
for  the  two  problems,  as  extracted  from  M/A  by  the  methods  of 
the  previous  section,  are  shown  in  Figures  3.1  and  3.2,  re- 
spectively. Since  M/A  does  not  calculate  ozone  profiles,  an 

[81] 

analytic  ozone  profile  taken  from  Green, 
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Figure  3.2>  — Test  Problem  Water  Vapor  Density  Profiles 
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(Z-Z  )/h 
We  ^ 

" ~1  (Z-Z  )/h-]i 

^ h [l  + e J 

was  used  with  parameters  W = 0.218  atm-cm  (total  ozone) , 
h = 4.5  I<m  (ozone  layer  half-width),  and  Z^  = 23.25  km  (ozone 
maximum).  This  profile  is  shown  in  Figure  3.3. 

The  cosine  of  the  solar  zenith  angle,  ground  tempera- 
ture, Tg  , and  surface  albedo  for  each  problem  were  also  taken 

directly  from  the  M/A  code.  For  Problem  1 they  are 

Vi  = 0.84641818 

T = 343.26°K 

g 

0 \)  < 2500  cm"' 

■ ) 

^0.2  V > 2500  cm"' 

and  for  Prcblem  2, 

0.454845 
304.19°K 

j 0 V £ 250 J cm" ' 

^0.09  V > 2500  cm"' 

Note  from  Figure  3.1  that  in  both  problems  there  is  a tempera- 
ture slip  at  the  ground,  of  magnitude  31.5°K  for  Problem  1, 
and  of  magnitude  15 °K  for  Problem  2. 

Both  problems  were  run  with  39  zones,  extending  from 
the  surface  to  2 mb.  A total  of  116  frequency  groups  wer< 
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used,  32  in  the  IR  (60-2500  cm"*)  29  in  the  near  IR  (2500- 
11,000  car  M and  55  in  the  "visible"  (11,000-50,000  cm  M • A 
variable  number  of  Gaussian  angles  was  used:  6 between  60  and 

800  cm"*,  10  between  800  and  1400  cm"*,  8 between  1400  and 
2500  cm”*,  and  10  between  2500  and  50,000  cm  *.  Fewer  angles 
were  used  in  the  strong  water  vapor  and  CO2  absorption  bands, 
based  on  sample  calculations  showing  a high  degree  of  hemi- 
spherical isotropy  in  the  intensity  in  these  bands.  Test  cal- 
culations were  run  with  more  angles  in  selected  frequency 
groups  to  determine  if  the  number  of  angles  being  used  was 
sufficient;  in  no  case  did  the  fluxes  vary  by  more  than  a 
percent  upon  using  more  angles. 

Test  calculations  were  also  performed  to  ascertain 
if  the  fluxes  were  sensitive  to  variations  in  the  exponen- 
tial fitting  process.  As  we  pointed  out  in  Section  2.5, 
the  exponents  and  coefficients  resulting  from  this  process 
are  highly  sensitive  to  perturbations  in  the  transmission 
function  data,  but  the  conjecture  was  made  that  the  fluxes 
and  intensities  would  be  insensitive  to  variations  in  the  fit. 

In  the  test  calculations,  various  small  perturbations  in  the 
transmission  function  data,  which  produced  large  perturbations 
in  the  exponential  fits,  nevertheless  produced  only  small 
perturbations  (always  less  than  1 percent)  in  the  fluxes. 

Hence,  our  previous  conjecture  has  been  borne  out  remarkably 

well. 

The  ATRAD  frequency-integrated  radiative  flux  F(z) 
for  Problem  1 is  given  in  Figure  3.4  as  a function  of  pressure 
p(z) . (NOTE:  positive  fluxes  are  downward  directed,  negative 

fluxes  are  upward  directed.)  The  Problem  1 atmosphere  is  very 
dry,  which  is  manifested  in  the  near  constancy  of  F throughout 
most  of  the  atmosphere.  The  large  change  in  F near  the  sur- 
face can  be  explained  in  terms  of  the  large  temperature  slip 
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Figure  3.4  - Total  (frequency-integrated)  radiative 
flux  profile  from  ATRAD,  Problem  1 
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(31.5°)  at  the  surface;  this  slip  produces  a large  upward 
blackbody  flux  which  decreases  the  downward  flux  and  which 
penetrates  several  hundred  meters  upward  before  absorption 
and  re-emission  by  water  vapor  can  produce  a cancelling  down- 
ward flux. 


The  Problem  1 flux  spectrum  F^(z)  for  z - 0 (top 
of  the  atmosphere)  and  z = z^  (surface)  is  plotted  as  a 
function  of  wavenumber  v in  Figure  3.5.  Actually,  because 
of  the  logarithmic  wavenumber  scale,  vF^  rather  than  F^  is 
plotted;  the  relationship 


/' 


vF  d(ln  v)  = /F  dv 


then  insures  that  the  area  between  the  curve  and  the  horizontal 
axis  in  any  wavenumber  interval  is  really  the  integrated  flux. 
Prominent  absorption  features  due  to  CO2  around  700  cm  ' (15y 

band),  O3  around  1000  cm"*  (9.6y  band),  H2O  in  the  near  IR 
(5  strong  bands  between  2500  and  12,000  cm  *) , and  0^  in  the 
ultraviolet  (35,000-50,000  cm"*)  are  all  visible  in  this 
spectrum. 

Comparisons  between  ATRAD's  and  Mintz-Arakawa  s radia~ 
tive  fluxes  and  "heating  rates"  are  presented  in  Table  3.1  for 
Problem  1.  The  "heating  rate"  table  simply  contains  flux 
differences  from  the  flux  table  above  it.  The  error  in  the 
M/A  net  heating  rate  is  quite  large  for  this  problem  - it  is 
too  small  by  about  a factor  of  2 for  the  lower  level,  and  too 
large  by  a factor  of  4 and  of  the  wrong  sign  for  the  upper 
level.  The  M/A  flux  into  the  ground  is  too  large  by  almost 
100  watts/m^ , or  about  31  percent.  The  M/A  fluxes  are  in 
general  considerably  more  accurate  than  the  M/A  heating  rates. 
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This  illustrates  the  fact  that  crude  models  to  obtain  radia- 
tive fluxes  perform  reasonably  well  in  practice  but  flux 
differences  obtained  from  such  models  can  be  grossly  in  error. 
It  is  necessary  to  have  quite  an  accurate  flux  computation  if 
flux  differences  are  desired,  for  in  the  process  of  differenc- 
ing one  or  even  two  significant  figures  of  accuracy  may  be 
lost. 

The  ATRAD  frequency-integrated  flux  F(z)  for  Problem  2 
is  given  in  Figure  3.6.  The  Problem  2 atmosphere  has  con- 
siderably more  water  vapor  at  all  levels  than  the  Problem  1 
atmosphere  — a factor  of  12  more  in  total  water  content  and  a 
factor  of  19  more  at  the  surface  — which  explains  the  slow  de- 
crease of  F with  altitude  (the  higher  in  the  atmosphere  one 
is,  the  stronger  the  upward  flux  due  to  water  vapor  emission 
becomes,  which  cancels  more  of  the  downward  solar  flux).  As 
in  Problem  1,  the  sharp  flux  change  near  the  surface  is  as- 
cribed to  the  temperature  slip  (15®) . The  change  is  much 
more  abrupt  in  Figure  3.6  than  in  Figure  3.4  because  the  air 
near  the  surface  is  so  much  wetter,  causing  the  water  vapor 
emission  to  cancel  the  temperature  slip  flux  a very  short 
distance  above  the  surface. 

The  Problem  2 flux  spectrum  F^(z)  is  given  in  Figure 
3.7  for  z = 0 and  z = z^.  For  the  same  reasons  as  in  Figure 
3.5#  vF^  rather  than  F^  is  plotted.  Note  that  the  near  IR 
H2O  absorption  features  are  enhanced,  and  that  the  CO2  15u 
and  Oj  9.6u  features  are  more  washed  out,  with  respect  to 
Problem  1.  This  is  because  of  the  increased  amount  of  water 
vapor  in  Problem  2. 

The  comparison  of  ATRAD  with  M/A  for  Problem  2 is  pre- 
sented in  Table  3.2.  The  M/A  and  ATRAD  fluxes  agree  much  more 
closely  for  Problem  2 than  they  did  for  Problem  1.  The  M/A 


3-15 


Pressure  (mb) 


vF  (cm".*  X watts/m 


Figure  3.7  _ Radiative  flux  spectrum  at  the  surface  and  at  the  top 
of  the  atmosphere.  Problem  2 
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net  flux  errors  are  1.6  percent  for  o = 0,  2.3  percent  for 
o = 0.5,  and  4 percent  for  a - 1.0.  Again,  hov/ever,  because 
of  loss  of  significance  in  subtraction,  the  heating  rate 
errors  are  substantially  larger.  The  M/A  heating  rate  in  the 
lower  level  is  too  large  by  more  than  a factor  of  2,  and  in 
the  upper  level  is  in  error  by  7.7  percent.  The  M/A  surface 
flux  (a  = 1) , which  determines  the  heating  of  the  ground,  is 
in  error  by  4 percent.  Thus,  the  M/A  radiation  model  makes 
very  good  predictions  of  fluxes  and  heating  of  the  lov;er 
level.  It  is  not  surprising  that  the  M/A  radiation  model 
performs  substantially  better  on  Problem  2 than  it  did  on 
Problem  1,  for  it  devotes  a great  deal  of  attention  to  water 
vapor  absorption,  and  Problem  2 is  much  more  dominated  by 
water  vapor  absorption  than  Problem  1. 

It  is  clear  that  the  M/A  radiation  model  incurs  sub- 
stantial error  in  the  case  of  dry  atmospheres.  However,  no 
definitive  conclusions  with  respect  to  the  usefulness  of  the 
M/A  radiation  model  for  global  climate  calculations  can  be 
drawn  from  just  these  two  problems.  Many  more  test  calcula- 
tions, including  cases  with  clouds,  cases  with  realistic 
aerosol  distributions,  cases  with  high  albedo  (such  as  in  the 
Arctic) , and  nighttime  cases,  need  to  be  performed  in  order 
to  obtain  a more  complete  picture. 

Even  when  this  picture  is  obtained,  other  factors 
need  to  be  considered.  For  example,  the  magnitude  of  the 
error  in  the  total  heating  rate,  including  latent  heating, 
adiabatic  heating,  etc.,  produced  by  errors  in  the  radiative 
contribution  must  be  ascertained.  Also,  it  might  be  the  case 
that  a radiation  model  whose  averaged  (zonally,  diurnally, 
etc.)  predictions  are  correct  is  sufficient  for  climate  dynam- 
ics calculations.  In  that  case,  it  would  be  the  average  of 
many  ATRAD  calculations  which  would  have  to  be  compared  with 
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the  average  of  many  M/A  radiation  calculations.  The  questions 
raised  in  this  paragraph  will  have  to  be  addressed  before  the 
amount  of  improvement  necessary  in  the  M/A  radiation  model 
can  be  ascertained. 

3.3  ARCTIC  STRATUS  PROBLEM 

The  choice  of  an  Arctic  stratus  problem  for  our  first 
cloudy  calculation  with  ATRAD  may  at  first  sight  seem  somewhat 
strange,  but  there  are  several  good  reasons  for  it.  First, 
we  were  involved  in  some  radiation  studies  for  AIDJEX  (Arctic 
Ice  Dynamics  Joint  Experiment)  upon  which  this  problem  bears 
directly.  Second,  Arctic  problems  come  very  close  to  satis- 
fying the  ATRAD  assumption  (see  Section  2.1)  of  horizontal 
homogeneity,  both  with  respect  to  the  surface  and  with  respect 
to  the  stratus  cloud  deck.  Obviously,  this  assumption  needs 
to  be  examined,  but  not  before  ATRAD  is  tested  in  situations 
which  closely  approximate  horizontal  homogeneity.  Third,  mea- 
surements were  made  on  three  different  days  with  widely  vary- 
ing albedos,  but  similar  cloud  conditions  so  that  a stringent 
test  of  ATRAD 's  response  to  albedo  variation  alone  was  possible. 
And,  localized  albedo  measurements  in  the  Arctic  can,  with 
greater  assurance  than  in  lower  latitudes,  be  taken  as  repre- 
sentative of  larger  areas  because  of  the  lack  of  substantial 
variations  in  the  terrain.  Finally,  cloud  size  distribution 
measurements  were  available,  which  is  unfortunately  not  always 
the  case  when  radiation  measurements  are  made. 

3.3.1  Comparison  of  ATRAD  with  Experiment 

The  measurements  to  which  we  will  compare  ATRAD  are 

[821 

reported  in  Weller,  et.  al.,‘  ^ and  were  made  near  Pt.  Barrow, 

Alaska,  during  the  month  of  June  1971  when  the  snow  was  melt- 
ing. The  data  consist  of  ground-level  solar  down-fluxes 
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(0.3  - 2.6y),  albedos,  and  solar  elevations  for  three  different 
days,  and  are  shown  in  Table  3.3.  Also  shown  are  the  flux 
measurements  "interpolated  to  30“"  according  to  Weller,  et.  al., 
and  the  corresponding  ATRAD  predictions  of  both  sets  of  fluxes. 


TABLE  3.3 

Comparison  of  flux  measurements  of  Weller,  et.  al.,^®^^  with 
flux  predictions  of  ATRAD.  (ATRAD  fluxes  have  been  rounded 
to  two  significant  figures  since  measurements  are  only  quoted 
to  this  accuracy)  . The  double— headed  arrow  indicates  the 
two  flux  values  whose  agreement  was  forced  by  adjusting  the 
ATRAD  cloud  droplet  number  density. 


Pt.  Barrow,  Alaska 

Date  (1971) 

June  2 

June  10 

June  29 

Albedo 

79% 

58% 

20% 

Solar  elevation 

31° 

mH 

21° 

Measured  surface 
down-flxix* 

0.56 

0.38 

0.26 

ATRAD  surface 
down  flux* 

0.60 

0.38 

0.26 

Measured  surface 
down-flux, * 
interpolated 
to  30°  solar 
elevation 

0.54 

0.47 

0.37 

ATRAD  surface 
down-flux, * 
30°  solar 
elevation 

0.58 

~ 

0.51 

0.43 

*A11  fluxes  in  cal/cm^/min 
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Before  discussing  Table  3.3,  let  us  specify  the  atmos- 
pheric structure  used  in  ATRAD  for  the  Arctic  stratus  problem. 
This  structure,  consisting  of  altitude,  pressure,  temperature, 
water  vapor  density,  ozone  density,  and  aerosol  (cloud)  number 
density,  as  taken  from  an  actual  ATRAD  run,  is  shown  in  Table 
3.4.  The  offset  of  levels  20  and  21  indicates  that  the  cloud 
layer  is  between  those  two  levels.  The  profiles  have  been 
taken  from  various  sources,  since  no  radiosonde  data  was  given 
in  Reference  [82]  . The  temperature  and  humidity  profiles  were 

taken  from  the  supplemental  Standard  Atmosphere  for  75°N  for 
[E3] 

July,  except  that  the  relative  humidity  values  within  and 
below  the  cloud  were  increased  to  100  percent.  Humidity  values 
above  10  km  were  obtained  from  the  stratospheric  polar  model 
of  Reference  [83]  , which  presumes  an  exponential  decrease  in 
mixing  ratio  by  a factor  of  10  between  10  km  and  12  km,  a con- 
stant mixing  ratio  between  12  km  and  17  km,  and  an  exponential 
increase  in  mixing  ratio  by  a factor  of  30  between  17  km  and  30 
km.  The  total  water  vapor  amount  is  1.65  g/cm*  of  which  0.33 
g/cm*  is  between  the  ground  and  the  cloud  base  and  0.30  g/cm^ 
is  within  the  cloud.  The  ozone  profile  usi2d  was  from  the  sub- 
arctic summer  atmosphere  of  McClatchey,  et.  al.,^^^^  resulting 
in  a total  vertical  ozone  amount  of  0.33  atm-cm.  The  cloud  is 
located  between  ^ km  and  1 km,  which  is  typical  for  Arctic  sum- 
mer stratus  according  to  Huschke.  ' The  cloud  size  distribu- 
tion is  taken  from  Weller,  et.  al.,  and  was  measured  near  Pt. 
Barrow  in  September  1971.  We  are  assured  that  these  measure- 
ments are  typical  for  summer  stratus,  however,  so  there  is  no 
difficulty  in  applying  them  to  our  June  situations.  The  size 
distribution  is  reproduced  in  Figure  3.8;  the  histogram  data, 
which  is  more  fundamental,  was  used  in  preference  to  the  smooth 
fit. 

The  ATRAD  calculation  used  22  levels  and  75  spectral 
intervals.  The  spectral  interval  structure  was  (in  cm~^) : 
3600(240)4800 (320)8000(500)32000(1000) 35000(1500)48500.  This 
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Figure  3. 8 - Number  distribution  of  Arctic  stratus  clouds  at  Barrow 
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structure  covers  the  instrumental  response  range  of  0.2  - 2.6p. 
Twelve  Gaussian  angles  were  used  in  each  spectral  interval. 

Since  the  atmospheric  structure  chosen  for  ATRAD  is, 
to  say  the  least,  somewhat  eclectic,  some  comments  need  to 
be  made  about  the  sensitivity  of  the  computed  fluxes  to  the 
input  data.  The  surface  fluxes  are  almost  totally  insensitive 
to  the  cloud  height  and  the  ozone  profile,  within  physically 
realistic  limits.  (The  measurements  of  Weller,  et. , al. , bear 
this  out  from  the  experimental  side,  as  far  as  the  insensi- 
tivity to  stratus  ceiling  height  goes)  . Sensitivity  to  the 
size  distribution  has  not  been  investigated;  however,  based 
on  the  success  of  calculations  which  go  so  far  as  to  approxi- 
mate Mie  scattering  by  isotropic  scattering,  and  on  the  success 
of  simple  models  such  as  we  shall  discuss  shortly,  one  gains 
the  impression  that  it  is  only  a few  gross  cloud  parameters 
(droplet  density,  mean  droplet  radius,  etc.)  which  are  really 
important.  It  might  even  reasonably  be  hypothesized  that  any 
size  distribution  which  leads  to  the  same  set  of  gross  cloud 
properties  will  also  lead  to  the  same  flux  predictions.  (ATRAD 
will  be  used  in  the  future  to  try  and  formulate  this  hypothesis 
more  rigorously)  . Since  we  are  dealing  only  with  the  solar 
spectrum  0.3  - 2.6y  (the  instrumental  response  range) , the 
temperature  profile  only  enters  the  calculation  through  the 
effective  absorber  amount  of  water  vapor  and  through  the 
Rayleigh  scattering  coefficient,  and  in  both  cases  the  derived 
quantities  are  insensitive  to  realistic  temperature  variations. 

The  input  data  to  which  the  calculation  ^ sensitive 
are  the  albedo,  solar  elevation,  moisture  profile,  droplet 
density,  and  cloud  thickness.  The  albedo  and  solar  elevation 
are  known  (see  Table  3.3).  The  moisture  profile,  while  not 
known,  cannot  differ  too  greatly  from  the  standard  one  that 
we  use,  for  Arctic  summer  conditions  are  considerably  less 
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variable  than  mid-latitude  and  tropical  regimes,  where  cumulus 
convection  and  frontal  systems  are  important.  And,  further- 
more, the  absorption  of  solar  radiation  by  water  vapor  is 
decidedly  a secondary  effect  in  this  problem  compared  to  the 
reflection  of  solar  radiation  at  the  cloud  top  and  bottom 
and  at  the  surface.  As  for  the  droplet  density  and  cloud 
thickness,  it  is  really  jus^  their  product,  which  is  propor- 
tional to  the  optical  thicknc'ss  of  the  cloud,  that  is  impor- 
tant. The  cloud  thickness  was  fixed  at  ^ km;  then  we  varied 
the  (unknown)  droplet  density  until  we  obtained  exact  agree- 
ment between  the  ATRAD  surface  down-flux  prediction  and  the 
measurement  in  Table  3.3  for  June  29.  The  droplet  density 
necessary  was  28.35  cm  Based  on  the  observation  in  Ref.  (82], 

that  the  cloud  was  similar  on  the  three  days,  the  same  droplet 
density  was  used  for  the  other  two  days. 

The  ATRAD  flux  prediction  for  June  10  agrees  exactly 
with  the  measurement  to  the  two  significant  figures  given  in 
the  measurement.  For  June  2,  the  ATRAD  flux  prediction  differs 
from  the  measurement  by  7%.  Considering  the  difficulty  of  the 
problem,  this  sort  of  agreement  is  remarkable  and  constitutes 
an  excellent  experimental  verification  of  ATRAD. 

The  source  of  the  7 percent  disagreement  on  June  2 
could  be  any  combination  of  several  factors.  One  factor 
of  course  is  errors  in  ATRAD  itself.  However,  the  disagreement 
could  also  be  eliminated  by:  (a)  changing  the  solar  elevation 

from  31®  to  29®,  (b)  decreasing  the  albedo  from  79  percent 
to  65  percent,  or  (c)  increasing  the  droplet  density  from 
28.35  cm  ^ to  43.5  cm  It  seems  highly  unlikely  that  the 
albedo  measurement  could  have  been  in  error  by  anything  like 
the  14  percent  necessary  to  completely  account  for  the  dis- 
crepancy; neither  does  it  seem  very  likely  that  the  droplet 
concentration  was  some  53  percent  higher  on  June  2 than  on 
the  other  two  days.  Therefore,  if  we  are  to  ascribe  the 
disagreement  to  measurement  errors  (other  than  errors  in  the 
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pyranometer  itself) , it  is  clear  that  (a)  is  the  most  likely 
candidate.  The  computed  surface  fluxes  are  particularly 
sensitive  to  the  solar  elevation  (and  surprisingly  insensitive 
to  the  albedo  and  the  droplet  concentration)  . The  measured 
values  of  solar  elevation  in  Reference  [82]  are  referred  to  as 
'mean  solar  angles',  implying  that  the  solar  elevation  varied 
during  the  course  of  the  measurement  and  that  subsequently 
some  sort  of  average  was  taken.  If  the  weighting  function 
used  in  this  average  was  improper  and  the  correct  average 
was,  say,  30®,  already  half  of  the  discrepancy  would  have 
been  accounted  for.  The  rest  of  the  discrepancy  could  then 
more  easily  be  accounted  for,  perhaps  by  a coiTibination  of  a 
too-high  measuied  albedo  and  a higher  cloud  optical  thickness 
on  June  2 than  on  June  10  and  June  29  (note  that  the  droplet 
cont  "ation  could  have  been  nearly  the  same,  but  with  a geo- 
mp  1 y thicker  cloud  on  June  2) . 

The  ATRAD  sensitivity  studies  of  the  last  paragraph 
have  an  important  bearing  on  any  radiation  measurements  which 
are  conducted  under  Arctic  summer  stratus.  They  indicate 
that  great  accuracy  is  not  required  in  either  the  albedo  or 
the  cloud  droplet  concentration  when  one  wishes  to  predict 
the  downward  radiation  flux  at  the  surface.  In  the  above 
examples,  an  18  percent  decrease  in  the  albedo  or,  alternatively, 
a 53  percent  increase  in  the  cloud  droplet  concentration,  were 
necessary  in  order  to  produce  a mere  7 percent  drop  in  the 
down-flux.  On  the  other  hand,  accurate  knowledge  of  the 
solar  elevation  0^  is  very  important.  This  is  because  the 
down-flux  above  the  cloud  varies  closely  as  sin  9^  , while, 
according  to  ATRAD  computations,  the  down-flux  at  the  surface 
has  the  even  stronger  dependence  (sin  6^)^'^^.  The  enhanced 
6g~dependence  of  the  surface  flux  is  due  partly  to  the  diffusion 
of  the  radiation  within  the  cloud  end  partly  to  the  strong 
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dependence  of  cloud-top  albedo  on  9 . At  any  rate,  because 

s 

of  this  sensitivity  to  6^  , care  must  be  exercised  in  the 
taking  and  processing  of  experimental  data  involving  even 
small  (1°  - 2°)  variations  in  sun  angle. 


3.3.2  Dependence  of  Surface  Flux  on  Solar  Elevation 


The  last  two  lines  of  Table  3.3  indicate  the  pitfalls 
of  extrapolating  data  taken  at  one  sun  angle  to  other  sun 
angles.  The  first  of  these  lines  contains  the  original  mea- 
surements 'interpolated  to  30°  solar  elevation and 
the  second  contains  the  ATRAD  predictions  for  9^  30°  , 

all  other  parameters  remaining  the  same.  While  we  cannot 
be  absolutely  certain,  it  appears  that  the  'interpolated' 
values  were  actually  extrapolated  according  to  an  assumed 
sin  0g  dependence,  as  witness  the  following  comparisons; 


0.56 

0754 


1.04 


0.38 

0.47 


0.81 


0.26 

0.37 


0.70 


sin 

31° 

sin 

30° 

sin 

24° 

sin 

30° 

sin 

21° 

sin 

"3^ 

1.03 

0.81 

0.72 


The  left-hand  ratios  are  of  measured  to  'interpolated'  fluxes, 
and  the  right-hand  ratios  are  of  the  sines  of  the  corresponding 
measured  6_'s  to  the  sine  of  30°.  It  is  obvious  from 
Table  3.3  that  ATRAD  does  not  predict  a sin  6 variation  of 

5 

surface  down-flux,  nor  should  such  a variation  really  be 
expected  except  for  a clear  sky.  Therefore,  an  extensive 
series  of  ATRAD  calculations  were  made  for  varying  albedo  and 
sun  angle  in  order  to  investigate  the  actual  functional  form 
of  the  extrapolation  to  other  sun  angles. 
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The  results  of  this  parameter  study  are  shown  in  Table 

3.5.  Both  the  full  down-flux  F+  and  the  purely  scattered 
s ^ 

part  F+  at  the  ground  are  included.  The  droplet  concentra- 


tion is  28.35  cm 
of  sin  6 , 


Suppose  that 


varies  as  some  power 


TABLE  3.5 

ATRAD  predicted  down-fluxes  at  the  grc  und  (F4')  for  the  full 
instrumental  range  3600  - 48500  cm“l  and  for^the  purely- 
scattering  spectral  range  11000  - 48500  cm~l,  as  a function 
of  albedo  and  solar  elevation  (6  ) 


0.58 


0.79 


F4 

^ -1 
3600  - 48500  cm  ^ 

(watts/m^) 

F4- 

5 -1 

11000  - 48500  cm 
(watts/m^ ) 

410.5 

299.6 

277.3 

203.3 

156.6 

13  5.5 

443.1 

324.1 

299.4 

220.0 

169.9 

125.5 

529.3 

389.5 

357.8 

264.3 

202.6 

150.5 

593.6 

438.6 

401.3 

297.7 

226.9 

169.1 
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F+  = F (sin  6^)  . (3. '' ) 

go  s 

Then  for  fixed  albedo,  each  pair  of  values  of  F+  can  be 

y 

used  to  eliminate  F and  solve  for  a.  It  is  found  that 

o 

a varies  in  the  narrow  range  1.50  - 1.58,  with  a tendency 

to  decrease  slightly  from  vl.56  for  6 e (30°, 40°)  to  vi.51 

for  e e(20°,30°).  The  best  over-all  value  of  a for 
s 

0ge(2O°,4O°)  is  1.53.  a is  furthermore  entirely  independent 
of  albedo.  Thus  the  entire  albedo  variation  resides  in  F^; 
however,  we  do  not  have  sufficient  data  to  empirically  fit 
this  variation.  Neither  do  we  know  the  dependence  of  either 
F^  or  a on  cloud  droplet  concentration.  Both  of  these 
dependencies  shall  be  investigated  in  the  future. 

If  a relation  similar  to  Equation  (3.1)  is  postulated 

for  just  the  purely  scattered  component  ®F+  of  Table  3.5 

exactly  similar  results  are  obtained.  The  -exponent  a varies 

between  1.49  and  1.56,  with  a best  over-all  value  of  1.51. 

It  decreases  slightly  from  1.54  in  the  30°  - 40°  range  of 

e to  1.49  in  the  20°  - 30°  range.  And  it  is  independent 
s 

of  albedo.  The  fact  that  the  exponent  is  almost  the  same  as 

for  the  full  solar  spectrum  is  surprising,  and  must  indicate 

that  the  scattering  is  the  relatively  dominant  influence  upon 

the  0 -variation,  with  absorption  playing  only  a minor  role, 
s 

3.3.3  Comparisons  of  ATRAD  with  Simpler  Models 

Weller,  et.  al.,  propose  a multiple  cloud-ground  reflection 
model  with  which  to  predict  their  experimental  results.  This 
model  is  practically  identical  to  that  used  for  cloudy  cases 
in  the  Mintz-Arakawa  and  for  that  reason  the  discussion 

below  takes  on  added  interest. 

The  model  of  Weller,  et.  al- , is  illustrated  in 
Figure  3.9,1.  The  cloud  albedo,  as  seen 


3-31 


SSS-R-75-2556 


^■c't 


e(l-a)3(l-a)*F^4’^ 


/ 


•''A 


aB(l-a)  U-a)F^*j. 


cloud 


Figure  3.9  -Multiple  cloud-ground  reflection 
model  of  Weller,  et.  al. I®2]  p A 
is  the  dovm-flux  at  the  cloud  top, 
a and  a are  related  to  the  cloud 
absorption  and  albedo,  respectively, 
and  fcs  is  the  surface  albedo. 


from  either  the  top  or  the  bottom,  is  equal  to  a(l-a) . The 
cloud  absorbs  a fraction  a of  the  down-flux  incident 

on  the  cloud  top.  3 is  the  surface  albedo.  No  absorption 
between  the  cloud  and  the  surface  is  assumed  and  therefore 
the  down-flux  F4-  at  the  surface  resulting  from  the  multiple 
reflection  is; 

00 

F+  = FA(l-a)  (1-a)  S la3(l-a)]’^ 

9 Ct 


(1-a) (1-g) 

1 - a3(l-a) 


(3.2) 


Unfortunately,  Reference  (82]  attempts  to  fit  this  model  to  the 
fluxes  which  were  incorrectly  extrapolated  to  30°,  arriving 
ti.ereby  at  parameter  values  a = 0.55  and  a = 0.07.  For 
completeness,  however,  the  predictions  with  these  parameter 
values  are  included  in  Table  3.6. 
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In  order  to  ascertain  the  predictive  capabilities  of 

Equation  (3.2)  as  a function  of  $ , for  fixed  sun  angle, 

the  parameters  a and  a were  chosen  to  give  the  ATRAD  values 

of  cloud  top  albedo  and  cloud  absorption  when  3=0, 

6 = 30°.  Referring  to  Figure  3.9,  it  may  be  seen  that: 

s 


F t (3  = 0) 
" F^4^(3  = 0) 


(1-a)  (1-a) 


F^(3  = 0) 

F^(3  = 0) 


(3.3) 


(3.4) 


The  ratios  on  the  right-hand  sides  were  taken  from  an  ATRAD 

calci’lation  with  3=0,  6_  = 30°,  yielding  for  a and  a , 

s 


and 


a = 0.4541  , 

a = 0.0274  . 


Using  these  values  of  a and  a , the  predicti«^ns  of  Equa- 
tion (3.2)  and  the  predictions  of  ATRAD  for  the  ratio 
F^/F^t  are  compared  in  Table  3.6  for  various  values  of  3. 
(The  other  columns  in  Table  3.6  shall  be  discussed  later). 

The  parameters  a = 0.55  , a = 0.07  , clearly  lead 
to  poor  predictions  vis  a vis  ATRAD.  The  parameters 
a = 0.4541,  a = 0.0274  of  course  lead  to  agreement  with 
ATRAD  for  3=0,  since  this  has  been  rigged,  but  as  3 
increases,  the  predictions  of  Equation  (3.2)  increasingly 
exceed  those  of  ATRAD.  Some  physical  effect  is  apparently 
being  ignored  in  the  simple  model.  Equation  (3.2). 
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TABLE  3.6 

Ratio  of  the  down-flux  at  the  ground  to  the  down-flux  at 
the  Arctic  stratus  cloud  top,  as  predicted  by  various 
models  for  solar  elevation,  30°,  for  several  values  of 
surface  albedo,  6 


6 

Eq.  (3.2) 
a=0.4541, 
a=0.0274 

Eq.  (3.2) 

a=0.55, 

a=0.07 

ATRAD 

Eq.  (3.12) 
parameters 
from  Eq. ( 3 . 10) 

Eq.  (3.17) 
parameters 
from  Eq,  (3.21) 

0 

0.5310 

0.4185 

0.5310 

0.5310 

0.5310 

0.20 

0.5824 

0.4662 

0.5731 

0.5817 

0.5731 

0.58 

0.7138 

0.5950 

0.6748 

0.7106 

0.6749 

0.79 

0.8155 

0.7023 

0.7484 

0.8098 

0.7484 

There  are  two  physical  effects  which  an  examination  of 
the  ATRAD  fluxes  show  to  be  marginally  important.  One  is  the 
absorption  of  solar  radiation  within  the  cloud— to-ground 
(c'+'g)  layer  (by  the  near  IR  bands  of  water  vapor)  and  the 
other  is  the  Rayleigh  back-scattering  from  the  c-*-g  layer. 
Therefore,  a more  complex  multiple  reflection  model  accounting 
for  c->g  absorption  and  back-scatter  has  been  derived.  (Note 
that  the  a and  a derived  from  Equations  (3.3)  and  (3.4) 
contain  these  effects  in  some  crude  fashion,  since  the  ATRAD 
fluxes  do;  however,  until  the  effects  are  made  explicit,  the 
extent  to  which  they  are  accounted  for  can  only  be  conjectured) . 
In  the  process,  the  notation  of  Figure  3.9  has  been  rejected 
in  favor  of  the  following; 

= cloud  albedo, 

a = cloud  fractional  absorption, 

= albedo  of  c->g  layer, 

and 

= fractional  absorption  in  c-»-g  layer. 
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These  quantities  are  related  to  those  of  Figure  3.9  by: 
= a(l-a) , 

and 


We  also  define: 


= transmission  of  cloud 
= (1-a^) (1-a^)  , 

and 

= transmission  of  c-^g  layer 
' (l-V  (l-a„)  . 

Consider  now  the  case  3=0,  illustrated  in  Figure  3.10. 
The  upward-directed  fluxes  at  the  cloud  bottom  are  then  entirely 
ascribable  to  Rayleigh  back-scattering,  since  the  surface 
contributes  nothing.  If  we  sum  the  infinite  geometric  series 
implied  by  the  "etc.”  in  Figure  3.10,  we  arrive  at  the  equations: 


_ Fi(3=0) 

^1  ■ F^\(B=0)  " ^c'^a'^c  ' 


F 1(B=0) 

^2  = F^^TPoT  " “c  ^c“r^c^  ' 


cuid 


_ F^^(3=0) 
^3  = F^+^(3=0) 


^4  - F^+^(3=0) 


r T 

c c 


^cVc 


(3.5) 


(3.6) 


(3.7) 


(3.8) 
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where 


and  "cb"  refers  to  the  cloud  base 


Vc 

t'v 

I 

\ 

i 

i 

i 


a a„T 
c R c 


c \<  c 
\ 

\ 


T o^aiT 
a c R c 


cloud 


surface 

n'/  n 


Figure  3.10  - 


Multiple  cloud-ground  reflection 
model  including  the  effects  of 
absorption  and  back-scattering 
from  the  cloud-to-ground  layer 
for  the  case  of  zero  surface 
albedo,  3=0.  The  incident 
flux  at  the  cloud  top  is  taken 
to  be  unity. 


If  the  ratios  Y2'  ^3'  ^4  regarded  as  known  from  an 

ATRAD  calculation  with  3=0,  then  Equations  (3.5)  to  (3.8) 
uniquely  determine  the  four  parameters  u , a , a.,,  and  a , 
For  the  0g  = 30“,  3=0  ATRAD  calculation  of  the  Arctic 

stratus  problem,  the  parameters  are  found  to  be; 
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and 


a. 


a 


w 


a. 


'U 

— = 0.005687 
^3 


0.006850  , 


= 0.43998  , 

Y3  = 0.04229  . 


(3.10) 


Both  the  Rayleigh  back-scatter  and  the  absorption  a^ 

are  small;  they  could  only  be  important  when  the  surface 
albedo  is  high  and  the  radiation  bounces  between  cloud  and 
ground  many  times. 

Now  we  extend  the  model  to  8 > 0.  'In  order  to  avoid 
confusion,  all  the  multiple  reflections  between  the  cloud  and 
the  c->g  layer  are  summed  up  by  use  of  the  factor  of 

Equation  (3.9).  All  the  multiple  reflections  between  the 
ground  and  the  c-*-g  layer  are  summed  up  by  use  of  the  factor 


With  these  simplifications,  the  model  is  presented  in  Figure 
3.11.  The  algorithm  used  in  its  construction  is  as  follows: 

(a)  a down-flxix  at  the  cloud  base  coming  from  the 

cloud  top  is  enhanced  by  the  factor  F^,  then 

a fraction  a of  this  is  reflected  upward 
R 

and  proceeds  without  further  reflection  (but 
attenuated  by  the  factor  T^)  to  the  cloud  top; 


3-37 


SSS-R-75-25i36 


(b)  a down-flux  at  the  ground  is  enhanced  by  the 
factor  Fg  , then  a fraction  3 of  this  is 
reflected  upward  and  proceeds  without  further 
reflection  (but  attenuated  by  the  factor  T ) 
to  the  cloud  base; 

(c)  an  up-flux  at  the  cloud  base  coming  from  the 

ground  is  enhanced  by  the  factor  F^  , then 

a fraction  a of  this  is  reflected  downward 
c 

and  proceeds  without  further  reflection  (but 
attenuated  by  the  factor  T^)  to  the  ground. 

Note  that  part  (a)  of  the  algorithm  is  only  exercised  once, 
since  we  are  neglecting  Rayleigh  back-scattering  of  the  up- 
fluxes  at  the  cloud  top.  This  neglect  is  justifiable  because, 
first,  the  effect  is  small,  and  second,  the  error  is  partially 
compensated  by  looking  only  at  the  ratios  of  fluxes  to  the 
down- flux  at  the  cloud  top  • 

Summing  the  infinite  geometric  series  implied  in 
Figure  3. 11, we  arrive  at  the  model; 


F+  F F T T 

g _ q c a c 

Ft  1 - a 3F  F T2 
ct  c g c a 


(3.12) 


F t 
ct 

F t 
Ct 


3F 

a_  + a„F  ^ ^ ^ 


R c c 


1 - a 3F  F T" 
c g c a 


(3.13) 


and 


Ft  FT 

cb  _ c c 

, F t 1 - a 3^'  F T2  ' 
ct  c g c a 


(3.14) 


^cb 


3F  r^T^T 

a F T + g...^  ^ ^ ^ 

R c c 1 - a 3F  F T2 
c g c a 


(3.15) 
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That  Equations  (3.12)  to  (3.15)  reduce  properly  in  various 
limiting  cases  may  be  verified.  Obviously  for  3=0,  they 
reduce  to  Equations  (3.5)  to  (3.8).  For  a = 0 , Equation 
(3.12)  reduces  to: 

F+  (1-a J (1-a^) (1-a  ) 

g w c c 

^ct  1 - (1-a  ) a^3 

w c 

which  is  similar  to  the  simple  model,  Equation  (3.2),  except 
that  the  transmission  of  the  c-^g  layer  is  (1-a^)  rathei  than 
unity. 

Using  Eq.  (3.12)  with  parameter  values  from  Eq.  (3.10) 
leads  to  the  column  which  is  so  labeled  in  Table  3.6.  The 
flux  ratios  predicted  by  Eq.  (3.12)  are  slightly  lower  for 
3 > 0 than  those  predicted  by  the  simpler  model  of  Eq.  (3.2), 
and  the  adjustments  are  in  the  right  direction  to  bring  the 
ratios  into  agreement  with  ATRAD.  H'^wever,  the  adjustments 
are  much  too  small.  Hence  we  must  look  further  for  the  physi- 
cal effect  which  has  been  ignored. 

A fact  which  is  often  overlooked  vis  a vis  the  parameter 
which  we  call  'albedo'  is  that  this  quanta cy  is  not  an  intrin- 
sic property,  but  depends  on  the  angular  distribution  of  the 
incident  radiation.  This  is  true  whether  we  are  speaking 

of  surfaces  or  cf  clouds,  but  for  clouds  the  angular  dependence 
is  particularly  strong.  For  example,  for  the  Arctic  stratus 
cloud,  ATRAD  predicts  a cloud-top  albedo  of  '\'52  percent 
when  6 =20°  and  of  ^^^44  percent  when  9 =30°.  Hence 

it  is  unreasonable  to  expect  the  top  and  bottom  of  a cloud 
to  have  the  same  albedo,  for  they  experience  quite  different 
fields  of  radiant  intensity.  The  intensity  at  the  cloud  top 
is  primarily  collimated,  while  that  at  the  cloud  bottom  is 
thoroughly  diffuse.  In  order  to  account  for  this  effect,  we 
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shall  define: 


a,  = diffuse  albedo  of  cloud 
d 


and  use  for  the  albedo  of  the  cloud  bottom.  Defining 

d 

the  new  symbols: 


t = transmission  of  cloud  from  bottom  to  top 
c 


(1-a^) (1-a^) 


r.  = 


(3,16) 


the  revised  model  becomes: 


F + 

_SL_ 


r r ,T  T 
g d a c 

1 - a,er  r,T! 

d g d a 


F t 
ct 


ot 


a + 
c 


6r  r^T^T  t 

. qdacc 

1 - a.er  r;T-^ 


‘d*"  g d a 


(3.17) 


(3.18) 


TdTc 


(3.19) 


er  r^T^T 

r.  ^ q d a c 

r .T  + T 2 — 

R d c 1 - a.Br 


r .T2 

g d a 


(3.20) 


With  the  addition  of  , Equations  (3.17)  to  (3.20) 

are  now  a five-parameter  model.  Therefore  the  four  values  of 
the  above  ratios  at  6 = 0 Y2"  Y3#  Y4;  see  Equations  (3.5) 

to  (3.8))  are  insufficient  to  determine  all  five  parameters. 
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One  finds,  however,  that  a„  and  a are  given  by  the  same 

expressions  as  in  Equation  (3.10)  and  therefore  they  have 

the  same  numerical  values  as  given  there  (for  6^  = 30°)  . 

This  is  as  it  should  be  since  the  model  alteration  concerns 

the  cloud  and  not  the  air  beneath  it.  To  solve  for  a^,  a^, 

and  a^  , we  shall  take,  as  an  additional  datum,  the  value 

of  F4-/F  +.  from  an  ATRAD  calculation  for  = 30°,  3 = 0.20, 

g ct  s 

F+ (3=0.20) 

^5  Ft  (3=0.20)  • 

c t 

It  is  then  possible  to  deduce  that: 


= 0.005687  , 


a = 0.006850  , 

w 


= 0.3706 

a = 0.4398 
c 

and 

a = 0.04226 
c 


(3.21) 


Using  this  parameter  set  in  Equation  (3.17),  the  last  column 
of  Table  3.7  is  generated.  The  agreement  with  ATRAD  is  exact 
for  the  two  albedos  0.58  and  0.79  for  which  the  result  is  not 
rigged.  Thus,  the  model  seems  very  promising. 

In  order  to  test  the  model  further,  ATRAD  calculations 

were  made  at  solar  elevations  6 of  20°  and  40°.  Clearly, 

s 

we  expect  to  vary  with  6^.  If  the  model  is  correct, 

however,  we  should  expect  to  remain  practically  unchanged. 

Therefore,  we  hold  fixed  at  0.3706  and  recalculate  only 

the  other  four  parameters,  forcing  them  to  agree  with  the  3 = 

0 ATRAD  calculation  as  before.  For  0^  = 20°,  the  parameters  are 
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Comparisons  of  ATRAD  and  Eq.  (3.17)  predictions  of  flux  ratio 

(down-flux  at  the  surface  to  down-flux  at  the  cloud  top) 
for  solar  elevations  20°  and  40°,  for  various  albedos. 


F'*'/F  i 
g'  ct 

0 = 20° 

CD 

II 

o 

o 

s 

s 

Equation  (3.17) 

Equation  (3.17) 

Albedo 

and 

ATRAD 

and 

ATRAD 

Equation  (3.22) 

Equation  (3.23) 

0 

0.4578 

0.4578 

0.5984 

0.5984 

0.20 

0.4942 

0.4945 

0.G459 

0.6454 

0.58 

0.5824 

0.5835 

0.7603 

0.7586 

0.79 

0.6460 

0.6482 

0.8429 

0.8402 

a_  = 


a = 


= 


R 


a = 


w 


0.5165 

0.04443 

0.005705 

0.005562 


and  for  0^  = 40°  they  are 


a = 


a_  = 


a.,  = 


a = 


w 


0.3702 

0.03888 

0.005587 

0.007819 
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Using  these  parameter  sets,  the  multiple-reflection  model  (Eq, 
3.17)  predictions  of  the  ratio  of  surface  down-flux  to  cloud- 
top  down- flux  are  compared  with  the  correspond!  'TRAD  pre- 
dictions in  Table  3.6.  The  agreement  is  excellent.  The  differ- 
ences are  considerably  less  than  1%.  Thus,  the  multiple- 
reflection  model  (Eq.  3.17)  incorporates  all  the  major  physical 
effe  Its  and  can  be  relied  upon  to  extrapolate  the  surface  flux 
to  any  albedo,  if.  its  parameters  are  adjusted  to  predict  the 
correct  fluxes  for  6=0. 

It  is  of  interest  to  study  the  variations  of  the  param- 
eters a , a , a„,  a with  solar  elevation,  since  ultimately 
o c R w 

it  is  desirable  to  parameterize  these  variations  and  free  the 
multiple-reflection  model  entirely  from  its  dependence  on 
ATRAD.  Insufficient  calculational  data  prevents  our  making 
any  definitive  conclusions  as  yet,  however,  the  following 
trends  are  evident  for  typical  Arctic  solar  elevations: 

(a)  ttn  is  practically  independent  of  0„; 

(b)  a increases  slowly  with  6 , approximately 

/ • O nO.54  ® 

as  (sin9g)  ; 

(c)  a decreases  markedly  as  6 increases, 

c s 

and  its  variation  is  not  even  close 

to  a power  of  sin6  or  cos 6 ; 

s s 

(d)  a decreases  slowly  as  6 increases, 

c s 0 

very  approximately  as  (cosS^)^’ 

The  behavior  (c)  of  the  cloud-top  albedo  is  very  similar  to 
the  behavior  of  the  albedo  of  natural  surfaces,  including 
snow,  ice,  and  water.  Therefore,  empirical  formulas  which 
have  been  found  useful  for  fitting  the  albedo  of  natural  sur- 
faces should  be  applicable  to  a also*  The  behavior  of  a , 

c w 

which  should  be  strictly  constant  if  the  under-cloud  radiation 
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field  is  strictly  diffuse,  indicates  that  some  deviation  from 
diffuseness  exists  under  the  cloud.  And,  finally,  the  de- 
crease of  cloud  absorption  as  the  sun  rises  higher  is  under- 
standable in  terms  of  the  relatively  shorter  average  path  tra- 
versed by  a photon  impinging  on  the  cloud  top  at  a large  angle. 

3.4  ATRAD  COMPARED  WITH  KATAYAMA  MODELS  FOR  CLEAR-SKY  CASES 

In  Section  3.2  we  presented  comparisons  between  the 
predictions  of  ATRAD  and  those  of  the  older  Katayama  radiation 
model  as  used  in  the  2-level  Mintz*-Arakawa  GCM  at  RAND.  Tv;o 
clear-sky  problems,  a wet  one  and  a dry  one,  were  discussed. 
Further  analysis  has  subsequently  been  performed  on  those  prob- 
lems, and,  in  addition,  comparisons  have  been  made  with  the 
newer  Katayama  radiation  model  used  in  the  3-level  Mintz- 
Arakawa  GCM  at  UCLA.  The  "tter  model  has  been  found  to  give 
drcimatically  better  results  in  the  IR  as  compared  to  the  older 
version.  Problems  remain,  however,  in  the  parameterization  of 
near  IR  solar  absorption,  and  in  the  choice  of  a pressure  scal- 
ing factor. 

The  atmospheric  structures  for  the  two  clear  sky  prob- 
lems are  shown  in  tabular  form  in  Tables  3.8a  and  3.8b,  which 
were  taken  from  actual  ATRAD  runs.  (It  would  facilitate  inter- 
comparison of  models  if  all  authors  would  present  their  model's 
atmospheric  profiles  in  tabular  as  well  as  graphical  form. ) 

The  temperature  and  humidity  profiles  from  actual  M/A  (Mintz- 
Arakawa)  2-level  grid  locations  in  Chad  (North  Africa)  and 
Bolivia  were  interpolated  to  the  39  ATRAD  levels  according  to 
formulas  given  by  Gates,  et.  al.,^®^^  as  discussed  at  some 
length  in  Section  3.2.  The  Chad  and  Bolivia  problems  are  op- 
posite extremes  in  terms  of  water  vapor  content;  the  total 
vertical  water  vapor  amount  for  Chad  is  0.204  g/c^^,  while 
for  Bolivia  it  is  2.530  g/cm^,  or  12.4  times  as  much.  The 
ozone  profile  is  the  same  for  both  problems. 


Atmospheric  profiles  of  altitude  (km),  pressure  (mb),  temperature  (°K) , water  vapor 
(g/m’),  ozone  density  (atm-cm/km) , and  aerosol  density  (=  0)  for  Chad  problem. 
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TABLE  3.8b 

Atmospheric  profiles  of  altitude  (km),  pressure  (mb),  temperature  (°K) , water  vapor 
(g/mM f ozone  density  ( atm- cm/km ) , and  aerosol  density  (=  0)  for  Bolivia  problem. 
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(z-z^)/h 

atm- cm/km 


with  W = 0.218  atm-cm,  z^  = 23.23  km,  and  h = 4 . 5 km  . 

Due  to  insufficient  resolution  in  the  ozone  layer,  ATRAD 
actually  computes  a total  ozone  amount  of  0.200  atm-cm  rather 
than  0.218  atm-cm.  The  remaining  model  parameters  for  Chad 
and  Bolivia,  consisting  of  sun  angle,  albedo,  and  surface 
temperature,  are  given  in  Table  3.9.  Note  the  discontinuities 
in  temperature  at  the  surface  in  both  problems.  The  a-levels 
0,  1/2,  and  1 correspond  to  levels  11,  23,  and  39  in  the  ATRAD 
atmosphere. 

The  120  spectral  intervals  used  by  ATRAD  for  the  Chad 
and  Bolivia  problems  are  60(b0)  600 (20)  800  (40)  1200 (80)  1600  (160) 
2400(240)  4800(320)  8000(500)11000  (120)  11120  (380) 11500  (500)  32000 
(1000)35000(1500)48500  cm”^.  Anywhere  from  6 to  12  Gaussian 
angles  are  used,  depending  on  the  degree  of  isotropy  of  the 
intensity  in  each  hemisphere  (upward  and  downward; . 

A point  which  v;as  not  noted  in  the  previous  comparisons 
between  ATRAD  and  the  older  Katayama  ('OK')  model  was  that  the 
two  models  use  different  solar  constants.  ATRAD  uses  the  more 

[151 

current  value  of  1.94  ly/min  while  the  older  Katayama  model 
uses  2.00  ly/min.  l.iis  has  obvious  consequences  for  the  fluxes. 


TABLE  3.9 

Solar  elevation  (0g),  albedo  (6),  and  surface  temperature 
(Tg)  for  Chad  and  Bolivia  problems.  (B=0  for  X^3y) . 


Chad 

Bolivia 


0 

s 

B 

— 

Tg(°K) 

57.824 

0.20 

343.26 

27.055 

0.09 

304.19 

3-48 


and,  in  fact,  correcting  the  'OK'  solar  fluxes  to  the  more 
modern  solar  constant  produces  much  better  flux  agreement  be- 
tween the  two  models  in  the  solar  spectrum.  The  revised  flux 
tables  are  presented  in  Table  3.10.  In  spite  of  better  agree- 
ment in  the  solar  spectrum  (differences  no  larger  than  2.6% 
at  any  level),  the  poor  agreement  in  the  IR  still  causes  the 
net  fluxes  to  differ  by  as  much  as  2 0%.  The  average  net  flux 
disagreement  between  the  two  models,  for  all  three  levels  and 
for  both  problems,  is  8-1/2%.  At  their  face  value,  errors  of 
even  as  much  as  20%  may  not  seem  disturbing.  However,  it  is 
well  to  remember  that  we  have  been  discussing  fluxes  and  not 
flux  differences.  The  picture  as  regards  flux  differences  is 
much  more  bleak,  as  we  shall  see  below. 

Some  comments  on  Table  3.10  are  needed  here,,  in  light 
of  the  fact  that  these  fluxes,  and  others  like  them,  are  used 
to  construct  the  flux  difference  tables  to  follow.  The  first 
fact  to  note  is  the  accuracy  to  which  the  fluxes  are  quoted. 

It  has  been  found  possible,  in  varying  the  spectral  interval 
structure  used  in  ATRAD  for  the  Chad  and  Bolivia  problems,  to 
produce  changes  in  the  fluxes  in  the  first  place  after  the 
decimal  point.  Such  changes  are  ascribed  primarily  to  the  ex- 
ponential fitting  algorithm,  and  to  the  way  the  intervals  are 
arranged  relative  to  the  absorption  bands.  Perhaps  when  some 
of  the  problems  are  ironed  out  of  the  fitting  algorithm  (cf. 
Section  2.5)  this  limitation  on  the  accuracy  of  ATRAD  flux 
predictions  will  become  less  important,  in  the  meantime,  the 
ATRAD  fluxes  quoted  in  Table 3. 10  are  estimated  to  have  an  un- 
certainty of  ± 0.5  watts/m^  due  essentially  to  the  lack  of 
infinitely  fine  spectral  resolution.  The  impact  of  other 
errors  (such  as  in  the  McClatchey  transmission  datJ^^^  used  by 
ATRAD)  upon  the  fluxes  has  not  been  estimated  as  yet. 


Conparison  between  ATRAD  (unparenthesized)  and  older  Kataya 
thesized)  flux  predictions,  in  watts /m^ , for  Chad  and  Boliv 
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The  dividing  point  between  the  IR  and  solar  spectra  is 
not  specified  in  the  Katayama  models,  and  so  the  value  2640 
cm“'  (3.79u)  was  selected  based  on  the  observation  that  in  the 
interval  2400-2640  cm"*  the  ATRAD  net  fluxes  were  primarily 
upward  while  in  2640  — 2880  cm  ' they  were  all  dov^nward.  Moving 
the  IR  dividing  point  to  2880  cm  * (3.47y)  would  cause  a uni- 

form decrease  in  the  magnitude  of  the  IR  fluxes  of  about  1.7 
w/m^  for  Chad  and  between  1.1  and  l.b  w/m^  for  Bolivia.  Moving 
it  to  2400  cm“  * (4.17vi)  would  cause  an  increase  in  the  magni- 

tude of  the  IR  fluxes  of  less  than  0.9  5 w/m^  for  both  Chad 
and  Bolivia.  Of  course,  there  would  be  compensatory  changes 
in  the  solar-flux  in  both  cases.  The  other  dividing  point, 
at  11,120  cm“*  (0.8993y),  is  not  exactly  equal  to  the  Katayama 
separation  point  of  0.9y.  Again,  the  change  in  visible  vs. 
near  IR  flux  resulting  from  this  discrepancy  is  less  than 
0.5  w/m^.  These  figures  are  quoted  to  show  that  the  flux 
breakdown  is  not  sensitive  to  the  choice  of  either  dividing 
point. 

Not  only  the  'OK'  value  for  the  solar  constant  (2.00 

ly/min) , but  also  the  fractional  breakdown  of  th ' solar  flux 

between  X > 0.9y  and  X < 0.9y,  can  be  called  into  question. 

Since  either  fraction  determines  the  other,  consider  only  the 

X < 0.9u  fraction,  f . For  the  older  and  the  newer  Katayama 

o 

models,  f = 0.651.  From  the  detailed  spectral  data  of 
o '■151 

Thekaekara,  f = 0.634.  But  since  the  Katayama  models 

o 

neglect  ozone,  and  since  the  stratospheric  ozone  layer  will 

absorb  all  solar  flux  in  the  wavelength  region  X < 0.3y, 

which  amounts  to  1.2%  of  the  solar  constant,  it  should  be 

reasonable  to  take  f = 0.634  - 0.012  = 0.622  (while  leaving 

o 

the  X > 0.9y  fraction  at  1 - 0.634  = 0.366).  Or,  if  ozone  were 
to  be  included  in  the  same  way  that  stratospheric  water  vapor 
is,  the  X < 0.9y  part  of  the  solar  flux  could  be  appropriately 
attenuated  between  the  top  of  the  atmosphere  and  a = 0 (200  mb) , 
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in  which  case  f = 0.634  would  be  appropriate.  However,  be- 
o 

cause  of  the  way  in  which  absorption  of  solar  flux  is  param- 
eterized in  the  Katayama  models,  the  atmospheric  heating  rates 
in  the  solar  spectrum  are  independent  of  both  f^  and  the  solax 
constant.  Hence,  the  only  impact  of  changing  f^  and/or  the 
solar  constant  is  that  the  flux  into  the  surface  is  changed. 

Let  us  derive  the  appropriate  formulas  for  this  surface  flux 
from  the  treatment  in  Gates,  et.al.,^^^^  for  clear-sky  situa- 
tions . 

We  begin  by  defining  the  effective  water  vapor  amount  u* 

, . [80] 

between  the  surface  and  any  level  n: 

* 1 , ^/P_\“  dp  (3.24) 

n 

q is  the  water  vapor  mixing  ratio,  a is  the  pressure  scaling 
factor,  g is  the  acceleration  due  to  gravity,  is  the 

surface  pressure,  and  u*  must  be  in  units  of  g/cm^.  ?ot  the 
'OK'  model,  a = 1.  n = ® will  refer  to  the  top  of  the  at- 
mosphere. Also  define 


gsca 


f 


o 


S 

o 


y 


o 


gabs 


(1  - 


S u 

O O 


(3.25) 


which  are  the  'scattered'  and  'absorbed'  parts  of  the  solar 
flux  in  terms  of  the  fraction  f^  discussed  in  the  last  para- 
graph, the  solar  constant  , and  the  cosine  of  the  solar 

zenith  angle  u . The  fractional  absorption  of  the  'absorbed' 
o 

part  due  to  an  amount  u*  of  water  vapor  is  parameterized  by 
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A(u*,y^) 


0.189  , .0.303 

Vm_) 

(1  - f ) S ° 

o o 


(3.26) 


where  S is  the  mean  solar  constant  in  ly/min.  The  'absorbed' 
o 

flux  at  any  level  n is  therefore 


abs  ^ s^^^Tl  - A(u*  - u*,y  )1  (3.27) 

n " L n o J 

The  flux  differences  between  o = 0 and  a = 1/2  and  between 
o = 1/2  and  a = 1 are  then,  respectively, 


A 


1 


cabs 

^0 


abs 

2 


= 0.189 


y (s  /s  ) 

o o'  o 


ru 


^ - u* 


,0.303 


(3.28) 


A 


3 


cabs 

^2 


abs 

4 


= 0.189  y^ 


(3.29) 

The  ratio  (S  /^  ) depends  only  on  the  Earth-sun  distance, 
o o 

Therefore,  A^^  and  A^  are  independent  of  f^  and  . 

The  scattered  flux  at  all  levels  is 
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(1-3) (1-a  ) 
^sca  _ ^sca  ' ' o 

^ “ l-3a 


(3. 301 


where 


a = 0.085  - 


/P  (mb)  ^ 

0.247 


(3.311 


and  where  3 is  the  surface  albedo.  The  net  solar  flux  into 
the  surface  is,  therefore, 


S.  = 


gabs  ^ gsca 


(1-S) (1-a.) 


S y (d-f  ) [1-A(u*,y  ) ] + f 
O O ) O oo'  O 


l-3a 


(3.32) 


The  surface  flux  is  proportional  to  . (Lest  the  importance 


of  using  the  best  possible  value  of  to  calculate  be 

underestimated,  suffice  it  to  say  that  the  SMIC  Report  185] 


associates  a 1%  decrease  in  with  a decrease  in  global 


average  surface  temperature  of  j,.5'’C.)  The  dependence  of 


on  f^  is  slightly  more  complicated.  The  rate  of  change  of 


S.  with  f_  is 


dS. 


(1-3) (l-a„) 


= S y 
df  o o 


l-3oi 


- 1 


For  3e(0,1)  and  a^E(0,l)  , it  can  be  shown  that 


dS, 


3?^ 


< 0 
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Therefore,  if  we  decreased  f^  , as  proposed  in  the  last  para- 
graph, the  surface  flux  would  increase.  For  both  the 

Chad  and  Bolivia  problems,  was  calculated  from  Eq.  (3.32) 

for  f^  = 0.625  and  f^  = 0.651  in  order  to  determine  the 


quantitative  change 

in  S4 

produced.  The  results  for  Chad 

,827. 1 

w/m^ 

f = 0.651 
0 

^829.3 

w/m* 

f = 0.625 
0 

and  for  Bolivia, 

,423. 1 

w/m^ 

f = 0.651 
0 

^4  = 

M25.6 

w/m^ 

f = 0.625 
0 

The  change  in  in  both  cases  is  about  2 watts/m^ , or  0.3% 

for  Chad  and  0.6%  for  Bolivia.  This  change  is  produced  by  a 
4%  change  in  . Hence,  it  must  be  concluded  that  the  re- 

sults of  the  Katayama  models  are  insensitive  to  f^  , except 
perhaps  for  extreme  values  of  albedo,  sun  angle,  etc. 

While  flux  values  have  an  intrinsic  interest,  it  is 
really  only  flux  differences  between  atmospheric  levels 
which  are  used  in  a general  circulation  model  (the  only  ex- 
ception being  that  the  net  flux  at  the  surface  is  used  to  com- 
pute surface  heating).  In  our  experience,  the  ta)cing  of  such 
differences  almost  always  leads  to  the  loss  of  one  (and  some- 
times two)  significant  digits,  even  when  the  levels  are  as 
v^idely  sepa:ated  as  in  the  M/A  2-level  GCM.  Thus  a tlux- 
p.rediction  model  which  is  accurate  to,  say,  1%  will  have  un- 
certainties in  its  flux  difference  p-^edictions  of  the  order 
of  10%.  Even  though  the  'OK'  model  produces  reasonable  flux 
predictions,  therefore,  one  cannot  for  that  reason  alone  expect 
good  flux  difference  predictions. 
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In  Tables  3.11  and  3.12,  the  flux  differences  and  sur- 
face fluxes  as  computed  by  /TRAD  and  by  the  newer  and  older 
Katayama  models  are  presented.  ATRAD  computations  were  made 
with  and  without  the  ozone  profile  included,  in  order  to  assess 
the  influence  of  ozone  on  heating  rates.  Before  discussing 
these  tables,  however,  some  remarks  need  to  be  made  on  the 
newer  Katayama  ('NK')  model. 

The  'NK'  model  is  described  in  an  as  yet  unpublished 
manuscript  written  by  Katayama  at  UCLA.  Parts  of  this  manu- 
script as  well  as  a FORTRAN  listing  for  the  model  were  kindly 
furnished  to  us  by  Mr.  Hans  Giroux.  The  FORTRAN  version  of 
the  model  is  used  in  UCLA's  3-level  version  of  the  M/A  GCM. 

Rather  than  re-coding  the  'NK*  model  for  2 levels,  which  is 
non-trivial,  we  ran  it  'as-is'  by  putting  one  level  between 
a = 0 and  a = 1/2  and  two  levels  (o  = 1/2  to  a = 3/4  , 
a = 3/4  to  a = 1)  between  a = 1/2  and  a = 1 . This  puts 
the  incre  ised  resolution  where  most  of  the  water  vapor  resides. 

The  extra  temperature  values  required  (as  compared  to  the 
2-level  model)  were  obtained  using  the  temperature  interpolation 
formulas  of  the  'OK'  model.^®^^ 

Two  options  were  included  for  the  'NK'  humidity  profile. 
The  first  was  to  obtain  the  three  mixing  ratio  values  (Ql,  Q3, 

Q5)  required  by  the  'NK'  model  from  the  interpolation  formula 
used  in  the  'OK'  model. The  second  was  to  insert  the 
effective  water  vapor  amounts  (u*)  into  the  'NK'  model  directly, 
using  the  'OK'  values  of  q and  a (a  = 1)  in  Eq.  (3.24). 

The  first  option  will  be  called  the  'q-option,'  the  second  the 
'u*-option.'  Results  from  both  options  are  given  in  Tables  3.11 
and  3.12.  Since  both  options  use  Eq.  (3.24)  to  compute  u*, 
there  are  only  two  reasons  for  their  absorber  amounts  to 
differ:  (1)  they  use  different  q-profiles;  or  (2)  they  use 

different  pressure  scaling  factors  a . In  fact,  the  differ- 
ences in  absorber  amount  between  the  options  are  almost 
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Flux  differences  and  surface  fluxes  (watts/m^ ) for  Chad  as  predicted  by  ATRAD  with 
and  without  an  ozone  profile,  and  by  the  newer  cind  older  Katayama  models  ('NK'  and 
'OK'  respectively) . The  unparenthesized  'NK'  values  are  for  the  q-option,  while 
the  parenthesized  ones  are  for  the  u*-option.  The  upper  level  is  c = 0 to  a = 
1/2,  the  lower  level  is  a = 1/2  to  a = 1 . 
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entirely  due  to  different  a's.  (For  'OK',  a = 1;  for  'NK' , 

a = 0.6;  and  for  ATR?^D,  the  McClatchey  scheme uses  a = 0.9, 

1/2 

but  with  p/p  replaced  by  p/p  (T  /T)  ' to  account  for  the 
temperature  variation  oi  line  half-width.)  The  mixing  ratio 
profiles  used  in  the  'OK'  and  'NK'  calculations  were  identical 
up  to  pressure  level  p^^  (a  = 3/4)  in  both  models,  and  almost 
all  the  effective  water  vapor  amount  lies  below  this  level. 

The  'NK'  model  contains  no  improvements  in  the  param- 
eterization of  the  solar  spectrum  when  clouds  are  absent. 

Hence  the  'OK'  solar  flux  differences  agree  with  the  'NK'  u*- 
option  flux  differences.  The  'NK'  model,  however,  contains  a 
considerably  more  sophisticated  treatment  of  the  IR,  and  its 
predictions  of  IR  flux  differences  agree  much  more  closely 
with  ATRAD  than  do  the  'OK'  predictions.  In  the  Chad  problem, 
the  net  (solar  + IR)  flux  differences  ind  net  surface  flux 
predicted  by  'NK'  agree  well  with  ATRAD  and  are  a vast  improve- 
ment over  'OK'.  In  the  Bolivia  problem,  the  'NK'  and  'OK' 
net  flux  difference  predictions  are  close  to  each  other  in 
the  lower  level  (a  = 1/2  to  o = 1)  and  both  differ  by  roughly 
a factor  of  2 from  ATRAD.  In  the  upper  Bolivia  layer,  a 
fortuitous  cancellation  of  errors  in  the  solar  and  IR  flux 
differences  causes  the  net  'OK'  flux  difference  to  be  very 
close  to  the  ATRAD  value (s);  the  'NK'  net  flux  difference  is 
too  large  by  roughly  a factor  of  4/3.  Both  the  'NK'  and  'OK' 
surface  flux  predictions  are  off  by  about  14%  for  Bolivia. 

Thus,  the  following  general  conclusions  emerge  from  the  com- 
parisons in  Tables  3.11  and  3.12; 

(1)  for  a dry  or  a wet  atmosphere,  the  IR 
treatment  in  the  newer  Katayama  model 
is  far  superior  to  that  in  the  older 
model ; 
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(2)  V.here  is  no  clear-cut  advantage  in  the 
IR  of  the  q-option  over  the  u*-option, 
that  is,  of  a pressure  scaling  factor 
of  0.6  over  one  of  1.0;  in  the  (dry) 
lower  level  of  Chad,  0.6  ir.  clearly 
preferable,  while  in  the  (wet)  lower 
level  of  Bolivia,  1.0  is  clearly  pre- 
ferable; in  the  (dry)  upper  levels  of 
both  problems,  neither  0.6  nor  1.0  has 
a clear  advantage; 

(3)  the  solar  absorption  treatment  common 
to  both  Katayama  models  causes  too 
little  flux  to  be  absorbed  in  dry 
levels  (upper  and  lower  Chad, 

upper  Bolivia)  and  too  much  in  wet 
levels  (lower  Bolivia); 

(4)  for  a dry  atmosphere,  the  newer  Kata- 
yama model  is  vastly  superior  to  the 
older  one  in  predicting  net  heating 
rates  and  surface  fluxes;  for  a wet 
atmosphere,  the  newer  Katayama  model 
may  actually  be  slightly  worse  than 
the  older  one  in  predicting  those 
same  quantities; 

(5)  ozone  has  a substantial  impact  on 

the  heating  rates  of  the  upper  levels, 
and  even  affects  the  heating  rate  of 
the  Bolivian  lower  level  non-negligibly ; 
the  effect  is  concentrated  almost  en- 
tirely in  the  IR  and  thus  is  due  to  the 
9.6y  band. 
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In  order  to  quantify  the  statements  above  referring  to  'dry' 
and  'wet,'  let  us  note  that  the  effective  water  vapor  amounts 
for  Chad  and  Bolivia  are  as  given  in  Table  3.13.  Thus,  'dry' 
means  having  effective  water  vapor  content  less  than  ^ 0.2 
g/cm^  and  'wet'  means  having  effective  water  vapor  content 
larger  than  'v-  1.7  g/cm^. 

TABLE  3.13 

Effective  water  vapor  amounts  in  g/cm^  in  upper  (o  = 0 to  a = 
1/2)  and  lower  (o  = 1/2  to  o = 1)  levels  for  Chad  and  Bolivia, 
for  pressure  scaling  factors  (a)  of  0.6  to  1.0  and  for  ATRAD. 


CHAD 

BOLIVIA 

Upper 

Lower 

Upper 

Lower 

'NK' 

a = 0 . 6 

0.01715 

0.15294 

0.10267 

2.00645 

'OK' 

a = 1.0 

0.01217 

0.13775 

0.07249 

1.73746 

ATRAD 
a = 0.9 

0.01341 

0.13538 

0.07983 

1.76058 

In  view  of  the  fact  that  the  IR  seems  to  be  in  fairly 
good  shape  for  clear-sky  situations,  if  the  nBwer  Katayama 
treatment  is  used,  the  main  difficulty  in  predicting  clear-sky 
heating  rates  centers  around  the  absorption  cf  solar  radiation 
in  the  atmosphere.  In  order  to  localize  the  source  of  the 
error,  we  first  checked  the  Katayama  model  assumption  that 
Rayleigh  scattering  could  be  ignored  for  X > 0.9y.  For  the 
spectral  interval  3360  - 11,000  cm  ^ for  the  Chad  problem,  we 
found  that  ignoring  Rayleigh  scattering  increased  all  the  net 
fluxes  by  0.5  watts/m"  and  did  not  affect  the  heating  rates 


3-61 


SSS-R-75-2556 


at  all.  This  is  truly  a negligible  effect,  although  a more 
thorough  study  might  show  a marginally  non-negligible  effect, 
perhaps  at  high  albedos  and/or  low  sun  angles.  Similarly,  an 
examination  of  ATRAD  runs  with  and  without  ozone  (Tables  3.11 
and  3.12).  showed  that  the  largest  solar  flux  difference  change 
as  a resu].t  of  neglecting  ozone  was  0.9  watts/m  , which  cannot 
begin  to  account  for  the  errors  in  solar  heating  rate.  The 
only  other  important  absorbers  in  the  solar  spectrum  beside 
ozone  are  water  vapor  and  CO2.  Thus,  the  Katayama  model  errors 
must  be  due  to  some  combination  of  the  following  causes: 

(1)  the  parameterization  of  the  absorption 
of  solar  down-flux  by  wate  ; vapor  [Eq. 

(3.26) ] ; 

(2)  the  neglect  of  the  absorption  of  the 
surface-reflected  solar  up-flux; 

(3)  the  neglect  of  CO2  absorption. 

A comparison  of  the  solar  up-  and  down-flux  differences  for 
Chad  and  Bolivia  is  presented  in  Table  3.14,  we  have  broK. 
these  results  down  further  into  X < 0.9y  and  X > 0.9y  values, 
but  this  breakdown  is  rather  artificial  since  as  we  have  shown 
in  Eqs.  (3.28)  and  (3,29)  the  Katayama  parameterization  is 
independent  of  any  partitioning  of  the  solar  spectrum.  The 
breakdown  is  shown  primarily  to  illustrate  the  error  incurred 
if  one  were  to  lump  all  absorption  into  X > 0.9y. 

The  salient  facts  which  emerge  from  Table  3.14  are: 

(1)  there  is  a fundamental  disagreement  between  ATRAD  and 
the  Katayama  parameterization  of  absorbed  solar  flux  [Eqs. 
(3.28)  and  (3.29)]  whether  or  not  we  include  the  absorbed 
near  IR  (X  > 0.9y)  up-flux  and  absorbed  'visible'  (X  < 0.9y) 
flux;  (2)  the  absorbed  near-IR  up-flux  is  at  most  an  3% 
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Flux  differences  (watts/m^)  in  the  near-IR  (2640  - 11,120  cm"  ) and  'visible 
(11,120  - 48,500  cm“^)  for  the  Chad  and  Bolivia  problems.  Near-IR  values  ar 
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effect  (in  the  Chad  lower  level)  but  may  be  decidedly  non-neg- 
ligible  for  dry  problems  with  albedos  substantially  larger  than 
the  0.20  of  Chad;  and  (3)  the  absorbed  'visible'  flux  is  at  most 
an  11%  effect  (in  the  Bolivian  lower  level)  and  is  thus  marginally 
negligible.  Clearly  there  are  fundamental  difficulties  with 
the  function  A(u*,y  ) of  Eg.  (3.26);  it  simply  causes  too 
much  absorption  in  wet  levels  and  not  enough  in  dry  ones.  The 
simplified  model  which  we  shall  present  below  offers  an  at- 
tractive alternative  in  the  search  for  a replacement  for  A. 

Also,  we  shrill  comment  on  the  possibility  of  simply  improving 
on  A itself,  by  changing  its  parameters. 

The  simple  model  to  which  we  refer  ignores  Rayleigh 
scattering  and  thermal  emission  in  the  spectral  region  2640  - 
11,120  cm“*.  ATRAD  calculations  have  demonstrated  that 
including  Rayleigh  scattering  in  this  region  produces  a very 
small  net  flux  decrease  at  all  levels.  The  Planclc  function 
effect  is  even  smaller,  even  for  the  high  temperatures  of 
Chad.  Hence,  the  radiative  transfer  equation  reduces  to 


ve[2640  - 11,120] 


where  is  the  absorption  coefficient.  If  the  boundary 


condition  at  the  top  cf  the  atmosphere  is 


z=0 


= S 6 (^-^  ) 
V o 


(for  solar  flux  S^)  then  the  down-flux  is  simply 


- r a' 

^o  Jo  ^ 


( z)  dz 
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If  an  albedo  6 and  diffuse  reflection  at  the  surface  z - 
are  assumed,  the  up- flux  is 

Ft  = 26/t(Zg)  ^ a;(z)dz  (3.34) 

where  is  the  exponential  integral  of  third  degree.  Equa- 

tions (3.33)  and  (3.34)  constitute  a simple  model  for  the 
near-IR  when  the  atmosphere  is  aerosol-free.  Calculations 
have  so  far  been  performed  only  with  Eq.  (3.33)  in  order  to 
compare  both  with  ATRAD  and  with  the  Katayama  model. 

The  calculations  for  Eq.  (3.33)  were  actually  performed 
with  the  frequency-averaged  form 

H-,0  CO^r 

= Vav  [w*(0-z)/yoj  (3.35) 


where 


S 


Av 


Av 


/, 


v+Av 


S dv 
V 


(as  estimated  by  cubic  interpolation  in  Thekaekara's  tables 
and  where  is  the  transmission  function  for  either  water 

vapor,  for  vertical  H2O  amount  u*(0->-z)  between  0 and  z,  or 
for  vertical  CO^  amovmt  w(0->-z),  taken  from  McClatchey,  et. 
al.^^^^  Using  an  interval  Av  = 20  cm  ^ , the  fluxes  from  Eq. 
(3.35)  were  summed  across  the  entire  spectral  region  2640  - 
11,120  cm*"^  to  yield  the  results  of  Table  3.15.  Calculations 
were  made  with  and  without  CO2.  Comparison  values  of  down- flux 
and  down-flux  difference  are  presented  as  predicted  by  ATRAD 
and  by  Katayama  (Eqs . ( 3 . 26)  and  (3. 29) ) . The  ATRAD  spectral 
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interval  structure  for  these  calculations  was  2640(240)  4800 
(320)  8000(500)  11,000(120)  11,120  cm"^,  so  that  the  ATRAD 
frequency  resolution  was  considerably  coarser  than  that  of 
Eq.  (3.35).  The  predictions  of  Eq.  (3.35)  confirm  the  ATRAD 
predictions  remarkably  well,  and  in  fact  Table  3.15  furnishes 
an  excellent  validation  of  the  exponential  fitting  idea  upon 
which  ATRAD  is  based.  The  residual  differences  between  ATRAD 
and  Eq.  (3.35)  must  be  partly  due  to  Rayleigh  scattering, 
partly  due  to  the  approximation  of  T^^  by  an  exponential 
sum  in  ATRAD,  and  partly  due  to  ATRAD 's  coarser  spectral  reso- 
lution. Removing  CO-  from  the  atmosphere  causes  the  near-IR 

^ 1 

down-flux  to  change  by  at  most  1 - l-j%,  and  the  largest 
tropospheric  change  in  down-flux  difference  is  0.8  watt/m^ . 

The  primary  impact  of  CO2  in  the  near-IR  is  to  change  the 
stratospheric  heating,  which  is  grossly  mis-estimated  if 
only  water  vapor  is  considered.  Note  that  the  inclusion  of 
CO2  in  the  model  of  Eq.  (3.35)  actually  decreases  whe  heating 
rate  due  to  near-IR  down-flux  in  the  upper  and  lower  Bolivian 
and  lower  Chad  levels.  This  points  up  the  unreliability  of 
simple  intuition  in  complex  radiative  transfer  problems.  The 
only  reliable  intuitive  deduction  is  that  the  fluxes  must  de- 
crease when  the  effects  of  CO2  are  added,  and  Table  3.15  of 

course  bears  this  out.  In  the  discussion  of  the  solar  absorp- 

tion function  A which  follows,  the  Eq.  (3.35)  predictions 
without  CO2  shall  be  used  as  br..;hmark  values. 

None  of  the  results  of  Table  3.15  contradict  the  ob- 
servation, made  in  connection  with  Table  3.14,  that  the 
function  A [Eq.  (3.26)]  is  deficient  in  both  wet  and  dry 
situations.  An  obvious  direction  one  might  take  in  order  to 
improve  the  parameter ‘ "ation  would  be  to  study  the  predictions 
of  the  simple  model  Eq.  (3.35)  and  derive  a new  empirical 

function  to  replace  A.  In  the  long  run,  th^s  will  be  the 

only  satisfactory  way  to  proceed.  Undoubtedly,  this  new 
function  should  be  of  a tabular  nature  in  order  to  r.liminate 
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the  computing  costs  associated  with  analytic  functions  like 
Eq.  (3.26). 

In  the  shorter  term,  hov/ever,  it  would  be  desirable  to 
re-pararujterize  A in  order  to  make  it  more  accurate.  There 
are  three  possible  ways  to  do  this:  (1)  change  the  definition 

of  u*;  (2)  change  the  coefficient  (0.187);  and  (3)  change  the 
exponent  (0.303).  As  far  as  changing  u*  goes,  although  there 
are  infinitude  of  exotic  formulas  one  might  consider,  we 
shall  restrict  ourselves  to  observing  what  change  in  the 
pressure  scaling  factor  a (see  Eg.  (3.24))  would  improve 
matters.  But  we  already  have  sufficient  information  at 
hand  in  Table  3.11  and  3.12  to  observe  how  changing  a from 
1.0  (parenthesized  'NK'  values)  to  0.6  (unparenthesized  'NK' 
values)  affects  the  absorbed  solar  flux.  The  conclusion  one 
IS  forced  to  draw  is  that  there  is  no  net  gain  with  a = 0.6  - 
the  results  in  both  the  Chad  and  Bolivia  upper  levels  are 
slightly  improved,  in  both  lower  levels  slightly  worse. 
Furthermore,  the  changes  in  Katayama's  model  solar  absorption 
resulting  from  varying  a are  much  too  small  to  offer  any  hope 
that  agreement  with  ATRAD  can  be  obtained  in  this  fashion. 

Next  we  consider  if  A can  be  improved  by  changing 
its  exponeiit.  To  this  end,  we  present  in  Table  3.16  cal- 
culations oi  the  upper  and  lower  level  flux  differences 
(A^  and  A3)  from  Eqs.  (3.28)  and  (3.29)  with  varying  expo- 
nent, The  values  of  u*  used  are  those  of  the  'OK'  model  in 
Table  3.13.-  For  Chad,  a value  of  the  exponent  between  0.25 
and  0.303  would  lead  to  agreement  of  A^  with  Eq.  (3.35), 
but  the  value  of  A3  would  be  made  worse  thereby.  On  the 
other  hand,  there  is  no  exponent  which  would  lead  to  agree- 
ment of  A3  with  Eq.  (3.35);  the  maximum  A3  attainable  by 
varying  the  exponent  (36.3  watts/m^ ) is  still  12%  below  the 
Eq.  (3.35)  value.  As  for  Bolivia,  there  is  no  exponent 
which  will  bring  into  agreement  with  Eq.  (3.35);  no 
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TABLE  3.16 

Predictions  of  Eqs.  (3.28)  and  (3.29)  for  solar  flux  differ- 
ences (watts/mM  across  upper  (A]^)  and  lower  (A3)  layers  as 
we  very  the  exponent  from  the  Katayama  value  of  0.303,  for 
Chad  and  Bolivia  problems.  The  predictions  of  Eq.  (3.35) 
are  included  for  comparison. 


CHAD 


BOLIVIA 


Exponent 

‘3 

Exprinent 

‘1 

"3 

0.10 

29.2 

21.3 

0.10 

23.8 

19.4 

0.15 

31.7 

27.6 

0.15 

28.4 

28.9 

0.20 

30.8 

31.8 

0.20 

30.4 

38.3 

0.25 

28.1 

34.4 

0.215 

30.6 

41.2 

0.303 

24.6 

35.9 

0.25 

30.7 

47.9 

0.325 

23.1 

36.2 

0.303 

29.9 

58.0 

0.35 

21.4 

36.3 

0.325 

29.4 

62.3 

0.375 

19.7 

36.3 

0.35 

28.7 

67.2 

0.40 

18.1 

36.1 

0.40 

27.0 

77.1 

0.45 

15.1 

35.4 

0.50 

23.4 

97.8 

0.50 

12.6 

34.3 

0.60 

19.9 

120.0 

0.60 

8.5 

31.4 

Equation 

(3.35) 

27.7 

41.6 

Equation 
(3. 35) 

36.1 

42.1 
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matter  how  we  vary  the  exponent,  remains  at  least  15% 

too  small.  However,  an  exponent  of  about  0.215  will  force 
into  agreement  with  Eq.  (3.35),  and  this  exponent  also 
leads  to  a slightly  better  value  of  A^  . But  an  exponent 
of  0.215  leads  to  definitely  worse  values  of  and  A^ 

for  Chad.  Hence,  there  seems  to  be  no  exponent  which  is 
consistently  better  than  0.303. 

Next,  we  ask  if  any  value  of  the  coefficient  of  A 
would  be  better  than  0.189.  But  the  error  we  are  dealing 
with  is  not  monosigned.  A larger  coefficient  is  needed  for 
levels  with  a small  cimount  of  water  vapor,  and  vice  versa. 

Finally,  there  remains  the  possibility  that  some 
judicious  variation  of  both  the  coefficient  and  exponent 
might  improve  A . The  exponent  first  must  be  changed  so 
that  all  the  errors  are  monosigned,  then  the  coefficient 
changed  to  bring  A into  agreement  with  Eq.  (3.35).  An 
examination  of  Table  3.16  reveals  that  A^  for  Chad  can 
never  be  larger  than  36.3  watts/m^  (and  needs  to  be  41.6 
watts/m^)  and  that  A^^  for  Bolivia  can  never  exceed  30.7 
watts /m^  (and  needs  to  be  36.1  watts/m^)  no  matter  how  we 
vary  the  exponent.  Thus,  the  exponent  can  be  varied  so 
that  all  values  of  A^  and  A^  are  less  than  Eq.  (3.35) 
predictions,  but  not  vice  versa.  In  order  to  make  A^ 
for  Bolivia  less  than  42.1  watts/m^ , the  exponent  must  be 
less  than  0.215.  But  for  exponents  < 0.215,  we  find 

— ^3  (toughly)  for  Chad  so  that  no  coefficient  adjustment 
could  ever  produce  agreement  with  Eq.  (3.35). 

Thus,  the  resolution  of  the  difficulty  really  awaits 
the  creation  of  the  new  tabular  function  of  which  we  spoke 
earlier. ' 

3*5  ATRAD  COMPARED  WITH  KATAYAMA  MODEL  FOR  ARCTIC  STRATUS 

PROBLEM 

The  Arctic  stratus  problem  was  discussed  in  detail 
in  Section  3.3.  We  continue  that  discussion  here  by  briefly 
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comparing  some  of  the  Katayama  model  cloudy-sky  approxima- 
tions against  ATRAD. 

3.5.1  Scattering  in  the  IR 

It  is  assumed  in  the  Katayama  model,  and  indeed  by 
practically  all  extant  radia^ton  models,  that  scattering  in 
the  IR  can  be  neglected.  However,  until  the  advent  of  ATRAD, 
no  model  was  capable  of  actually  checking  this  assumption, 
for  the  combined  line  absorption-scattering  problem  was 
viewed  as  insuperable.  As  a start  towards  examining  this 
assumption  quantitatively,  therefore,  two  ATRAD  calculations 
were  made  for  the  IR  (60  - 1920  cm“M  spectral  region  for 
the  Arctic  stratus  problem.  The  first  included  Mie  scatter- 
ing in  the  normal  fashion,  the  second  was  identical  to  the 
first  in  all  respects  save  that  the  Mie  scattering  coeffi- 
cient (but  not  the  Mie  absorption  coefficient)  was  set  to 
zero.  Selected  results  from  these  calculations  for  down-flux, 
up-flux,  net  flux,  and  heating  rate  (flux  difference)  are 
presented  in  Table  3.17.  The  impact  of  ignoring  scattering 
may  be  summarized  as  follows; 

(a)  the  down-flux  below  the  cloud  is  decreased 
by  1.5  watts /m^ , or  0.5%;  above  the  cloud, 
it  is  unchanged; 

(b)  the  up-flux  above  the  cloud  is  increased 
by  1-2  watts/m^ , or  0.5%;  below  the  cloud, 
it  is  unchanged; 

(c)  the  net  flux  becomes  more  negative  by  1.1- 
2.2  watts/m^; 

(d)  heating  rates  are  unaffected  except  in  the 
cloud  (0.7%)  and  immediately  above  the 
cloud  (3.3%) . 
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TABLE  3,17 

Comparison  between  ATRAD  predictions  for  down-flux  FI,  up-flux 
Ft,  net  flux  F,  and  net  flux  differences  AF  (all  in  watts/m^) 
as  a function  of  altitude  z (km)  and  pressure  p (mb)  for  60  - 
1920  cm“*  for  the  Arctic  stratus  problem.  Unparenthesized 
values  are  for  full  Mie  treatment  and  parenthesized  values  are 
for  Mie  scattering  zeroed  out.  (The  cloud  is  located  between 
0.5  and  1.0  km,  T = 273°K  except  where  noted.) 


z 

P 

FI 

FI 

F 

0.0 

1012.5 

327.5  (326.0) 

311.7  (311.7) 

15. f ( 14.2) 

0.5 

951.9 

323.9  (322.3) 

322.4  (322. ^i) 

1.5  ( -O.i) 

1.0 

895.0 

236.7  (236.5) 

322.0  (324.0) 

-85.3  ( -87.5) 

2.0 

790.2 

210.1  (210.1) 

316.4  (318.0) 

-106.4  (-107.9) 

5.0 

536.7 

127.5  (127.5) 

244.8  (246.0) 

-153.7  (-155.0) 

10.0 

262.6 

43.7  ( 43.7) 

244.8  (246.0) 

-201.1  (-202.3) 

50.0 

1.0 

1.2  ( 1.2) 

245.5  (246.6) 

-244.3  (-245.4) 

1.0 

2.0 

5.0 

10.0 


10.0 

50.0 


AF 

AF  (T  =300°K) 
9 

-14.3  (-14.3) 

65.2 

-86.8  (-87.  4) 

-31.  3 

-21.1  (-20.4) 

- 20.7 

-47.3  (-47.1) 

-47.1 

-47.4  (-47.3) 

-47.2 

-43.2  (-43.1) 

-43.0 

SSS-R-75-2556 


r 


Thus,  the  neglect  of  Mie  scattering  is  borne  out  ex- 
tremely well  for  water  clouds  at  least  as  optically  thick  as 
the  Arctic  stratus  cloud.  The  strong  absorption  due  to 
liquid  water  in  the  IR  is,  of  course,  responsible  for  the 
suppression  of  scattering  in  the  cloud,  and  thus  this  result 
might  not  apply  to  an  aerosol  material  with  transmission 
v^indows  in  the  IR.  It  may  also  be  inapplicable  to  optically 
thin  clouds,  particularly  cirrus. 

3.5.2  Sensitivity  to  Surface  Temperature 

Because  of  the  .arge  inaccuracies  involved  in  the 
computation  of  the  Mintz-Arakawa  surface  temperature  T^, 
it  is  of  interest  to  assess  the  sensitivity  of  radiative 
heating  and  cooling  rates  in  the  atmosphere  to  this  parameter. 
A normal  ATRAD  calculation  of  the  Arctic  stratus  problem  for 
60  - 1920  cm“\  with  T^  = 273°K  was,  therefore,  compared  with 
a similar  calculation  with  T^  = 300°K.  In  the  first  case, 
the  surface  is  5°K  colder  tnan  the  cloud-to-ground  layer 
and  3"K  colder  than  the  cloud;  in  the  second  case,  the  sur- 
face is  22‘’K  warmer  than  the  cloud-to-ground  layer  and 
'v»  24°  warmer  than  the  cloud.  The  changes  in  heating  rate 
between  the  two  cases  may  be  observed  in  Table  3.16.  The 
cloud-to-ground  layer,  which  was  cooling  for  T^  = 273°,  is 
heating  for  T = 300°,  and  at  almost  five  times  the  rate  at 
which  it  was  formerly  cooling.  The  cloud  is  cooling  in 
both  cases,  but  only  one-third  as  rapidly  for  T^  = 300°  as 

for  T = 273°.  The  cooling  of  al3  levels  above  the  cloud 
g 

is  practically  unchanged,  so  that  the  cloud  effectively 
shields  .these  upper  levels  from  the  effects  of  a surface 
temperature  change.  In  spite  of  this  shielding,  however, 

AF  for  0 to  5 Km  is  -169.5  watts/m^  for  T^  = 273°K  and 
-33.9  watts/m^  for  T^  = 300°K,  so  that  the  cooling  of  the 
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whole  lower  half  of  the  atmosphere  is  dramatically  reduced 
(by  a factor  of  5)  as  a result  of  a 10%  increase  in  surface 
temperature . 

These  results  establish  the  sensitivity  of  IR  heat- 
ing and  cooling  to  surface  temperature.  There  is  very  little 
point  in  having  a reasonably  sophisticated  IR  treatment  in 
the  Mintz-Arakawa  GCM  if  it  is  only  to  be  coupled  with  an 
inadequate  predictive  scheme  for  T^. 

3.5.3  Cloud  Albedo 

A quantity  which  is  required  in  all  the  Katayama 

parameter! zations  relating  to  clouds  is  the  cloud  albedo, 

a . In  the  'OK'  model  a is  always  either  0.6  or  0.7.  In 
c c ■' 

the  'NK'  model  a assumes  different  values  for  X > 0.9y 

c 

and  X < 0.9y,  and  ranges  al3  the  way  from  0.19  (for  high 

cloud)  to  0.76  (for  cumulonimbus  cloud).  Foir  low  cloud, 

into  which  class  the  Arctic  stratus  cloud  must  fall,  the  'NK' 

model  assumes  a =0.66  for  X < 0.9y  and  a = 0.50  for 

c c 

X > 0.9y. 

Section  3.3  contained  an  extensive  discussi'-'n  of 

cloud  albedo,  particularly  as  it  related  to  surface  fluxes. 

A model  was  developed  for  cloud  'albedo'  (^ct'^^ct^ 

would  actually  be  observed,  as  a function  of  surface  albedo 

6 . Only  for  6=0  do  we  recover  the  actual  cloud  albedo 

, which  is  an  intrinsic  property  of  the  cloud,  from  a 

measurement  of  Ft . /F+.  (ct  = cloud  top) . The  Katayama  model 

ct  ct 

does  not  include  this  dependence  on  6 , nor  indeed  does  it 
include  any  dependence  on  droplet  concentration,  cloud 
thickness,  and  solar  elevation  0^  . ATRAD  is  an  ideal 
vehicle  for  studying  these  various  dependencies  with  a view 
to  parameterizing  them.  As  an  example,  consider  the  results 
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of  Table  3.18,  in  which  the  variation  of  ^ 

limited  range  of  sun  angles  and  albedos  is  shown.  It  is  im- 
portant to  remember  that  all  of  these  albedo  values  are  for 
exactly  the  same  cloud.  The  hopelessness  of  approximating 
clcud  'albedo'  as  a single  number,  rather  than  a function 
of  several  parameters,  is  brought  home  with  particular  force 
here.  We  feel  that  the  lack  of  a parameterization  for  cloud 
albedo  in  the  Katayama  model  is  its  most  serious  deficiency 
for  cloudy-sky  situations.  Until  such  a parameterization  is 
obtained,  the  other  formulas  in  the  Katayama  model  which  use 
c.' oud  'albedo'  can  only  be  believed  insofar  as  they  have  a 
certain  climatological  correctness  built  into  them. 


TABLE  3.18 

The  ATRAD-predicted  ratio  of  up-flux  at  the  cloud-top  (F+^) 
to  down- flux  at  the  cloud  top  (Fjt)  the  ArctuC  stratus 

problem,  for  various  values  of  solar  elevation  and  sur- 
face albedo  B.  Unparenthesized  values  are  for  the  full 
spectral  interval  3600  - 48,500  cm~\  The  first  parenthe- 
sized value  refers  to  3600  — 11,000  cm  only,  the  second 
parenthesized  value  refers  to  11,000  — 48,500  cm  only. 


bX 

O 

O 

CN 

30° 

O 

o 

0. 

0.20 

0.58 

0.79 

0.52  (0.52,0.52) 
0.58  (0.57,0.58) 
0.72  (0.69  ,0.74) 
0.83  (0.78,0.85) 

0.44  (0.44,0.44) 
0.51  (0.49,0.52) 
0.68  (0.64,0.70) 
0.80  (0.74,0.83) 

0.37  (0.36,0.38) 
0.45  (0.43,0.46) 
0.64  (0.59,0.66) 
0.78  (0.71,0.81) 
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3.5.4  Albedo  of  a Cloudy  Atmosphere 

Presume  that  we  know  the  cloud  albedo 


Pi 

ct 


for  X < 0.9y.  Then  let  us  enquire  into  the  accuracy  of  the 
Katayama  approximation 

a . = 1 - d-R'  ) (1-a^)  (3.36) 

atm  o 

for  the  albedo  a ^ of  the  whole  atmosphere  in  the  presence 

atm 

of  a single  cloud  layer,  is  the  albedo  of  the  correspond- 

ing clear  atmosphere  due  to  Rayleigh  scattering,  and  is  ap- 
proximated as  in  Eq.  (3.31)  by  the  Katayama  model.  In  Table 
3.19  are  assembled  values  of  for  the  Arctic  stratus 

problem  as  predicted  by  ATRAD  and  by  Eq.  (3.36)  for  the 
spectral  region  11,000  - 48,500  cm" d The  large  error  in 
the  formula  (3.36)  is  clearly  apparent.  In  these  examples, 

Eq.  (3.36)  seriously  overestimates  the  cloudy  atmospheric 
albedo,  whether  one  compares  it  with  the  ATRAD  value  at  the 
top  of  the  atmosphere  (1  mb)  or  at  167  mb  (the  difference  in 
the  ATRAD  values  at  1 mb  and  167  mb  is  due  to  Rayleigh  scat- 
tering and  ozone  absorption  between  the  two  levels) , Our 
parameter  studies  have  not  been  extensive  enough,  however,  to 
say  with  certainty  whether  Eq.  (3.36)  always  overestimates 
the  cloudy  atmospheric  albedo.  One  interesting  observation, 
however,  is  that  the  error  cannot  be  laid  on  the  approxima- 
tion (3.31)  for  . According  to  the  unpublished  Katayama 

manuscript  describing  the  new  model,  i:he  correct  should 

be  even  larger  than  Eq.  (3.31)  predicts  for  6^  = 33.6°. 
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TABLE  3.19 

Comparisons  of  cloudy  atmospheric  albedo  cx^^m  calculated 

by  Eq.  (3.36)  (using  the  R'  values  shown)  and  by  ATRAD  at 
1 mb  and  167  mb,  for  11,000  - 48,500  cm"’  for  the  Arctic  stra- 
tus problem  with  surfaoo  albedo  6=0  and  various  values  of 
solar  elevation  0g. 


0 

s 

R' 

(ATRAD) 

°‘atm 

(Eq.  3.36) 

“atm  <1  "*> 
(ATRAD) 

a (167  mb) 

atm 

(ATRAD) 

O 

o 

CM 

0.5162 

0.612 

0.497 

0.534 

30° 

0.4438 

0.531 

0.430 

0.459 

O 

O 

0.3758 

0.457 

0.369 

0.391 

Thus,  at  least  for  the  6 = 30°  case  in  Table  3.19,  the 

prediction  of  ^‘atm  would  be  even  worse  if 

the  correct  a were  used, 
o 

It  need  hardly  be  emphasized  that  significant  errors 

in  a have  serious  consequences  for  the  prediction  of 

atm 

climate.  Even  simple  global  climatic  models,  such  as  that 
of  Budyko,^^^^  point  up  the  sensitivity  of  the  Earth's  surface 
temperature  to  fluctuations  in  the  atmospheric  albedo. 


1 


r 

I 


REFERENCES 


Howell,  H.  and  H.  Jacobowitz,  "Matrix  Method  Applied  to 
Multiple  Scattering  of  Polarized  Light,"  J.  Atmos.  ScJ_^ 
(1970) , Vol.  27,  p.  1195. 

Dave,  J.  "Subroutines  for  Computing  the  Parameters  of  the 
Electromagnetic  Radiation  Scattered  by  a Sphere,  IBM 
Program  360D  - 17.4.002  (December  1968). 

Kerker,  M. , The  Scattering  of  lAght  and  Other  Electro- 
magnetic Radiation,  Academic  Press,  New  York  (1969). 

Hodkinson,  F. , Proceedings  oi  Inter-disciplinary  Confer- 
ence on  Electromagnetic  Scattering,  Potsdam,  N.Y_^,  lyt^, 
Pergamon  Press,  London  (1963) . 

Ellison,  J.,  Proc.  Phys.  Soc.  B,  (1957),  Vol.  10,  p.  102. 

Holland,  A.  and  G.  Gagne,  "The  Scattering  of  Polarized 
Light  by  Polydisperse  Systems  of  Irregular  Particles, 
Appl.  Opt.  (1970),  Vol.  9,  p.  1113. 


Hall,  F., 


;.  (1968)  , Vol.  7,  p.  891. 


Houghton,  J.  and  G.  Hunt,  "The  Detection  _ 

From  Remote  Measurements  of  Their  Emission  in  the  Far 
frared,"  Qtrly.  J.  R.  Met.  Soc.  (1971),  Vol.  97,  p.  1. 

Jacobowitz,  H.,  "A  Method  for  Computing  the  transfer  of 
Solar  Radiation  Through  Clouds  of  Hexagonal  Ice  Crystals, 
JQSRT  (1971),  Vol.  11,  p.  691. 

Kondrat'yev,  K. , Radiation  in  the  Atmosphere,  Academic 
Press,  New  York  (1969) . 

Hoover,  G.,  et.  al.  , "Infrared 

Bands  of  Atmospheric  Gases,  Appl.  Opt.  (-96  ),  • , 

p.  481. 

Milne,  E.  A.,  Thermodynamics  of  the  Stars,  in  Select^ 
Papers  on  the  Transfer  of  Radiation,  D.  Menzel,  Editor, 
Dover,  New  York  (1966).  ~ 


SSS-R-75-2556 


13.  Siegel,  R.  and  J.  Howell,  "Thermal  Radiation  Heat  Trans- 
fer, Vol.  Ill,"  SP-164,  NASA  (1971). 

14.  Woolf,  H.,  "On  the  Computation  of  Solar  Elevation  Angles 
and  the  Determination  of  Sunrise  and  Sunset  Times," 
TM-X-1646,  NASA  (1968). 

15.  Thekaekara,  M.  , et.  al.,  "The  Solar  Constant  and  the  Solar 
Spectrum  Measured  From  a Research  Aircraft,"  Technical  Re- 
port TR-R-351 , NASA  (1970). 

16.  Siegel,  R.  and  J.  Howell,  "Thermal  Radiation  Heat  Transfer, 
Vol.  I,"  Special  Publication  SP-164,  NASA  (1968). 

17.  Cox,  C.  and  W.  Munk,  "Slopes  of  the  Sea  Surface  Deduced 
From  Photographs  of  Sun  Glitter,"  Bull.  Scripps  Inst. 
Oceanog. , Vol.  6 (1956),  p.  401. 

18.  Boileau,  A.  and  J.  Gordon,  "Atmospheric  Properties  and  Re- 
flectances of  Ocean  Water  and  Other  Surfaces  for  a Low  Sun, 
Appl.  Opt.  , Vol.  5 (1966),  p.  803. 

19.  McStravick,  D. , "The  Effect  of  Surface  Roughness  on  the 
Reflected  and  Emitted  Energy  From  a Rough  Surface,"  Ph.D. 
Thesis,  Rice  University,  Texas  (1972) . 

20.  Paltridge,  T. , J.  Geophys.  Res.,  Vol.  76  (1971),  p.  2857. 

21.  Earing,  D.  and  J.  Smith,  "Target  Signatures,  Analysis  Cen- 
ter: Data  Compilation,"  Willow  Run  Labs.,  Institute  of 

Science  and  Techn.  , Univ.  of  Mich.,  Ann  A’'  jr  (1966). 

22.  Kyle,  T.,  "Calculations  of  Atmospheric  Transmittance  From 
1.7  to  20  U,"  JQSRT,  Vol.  9 (1969),  p.  1477. 

23.  Goody,  R.  , Atmospheric  Radiation,  Oxford  Univ.  Press, 
London  (1964)  . 

24.  Wyatt,  P.  , et.  al.,  Appl . Opt. , Vol.  3 (1964),  p.  229. 

25.  Smith,  W.  , "A  Polynomial  Representation  of  CO2  and  H2O 
Vapor  Transmission,"  ESSA  Tech.  Report  NESC-47,  National 
Environmental  Satellite  Service,  Washington,  D.C.  (1969). 

26.  McClatchey,  R. , et.  al. , "Optical  Properties  of  the  Atmos- 
phere," AFCRL-70-0527  , Air  Force  Cambridge  Research  Labs., 
Cambridge,  Mass.  (1970). 


R-2 


■tfMMliiillifeiWIli 


■H 


ik 


SSS-R-75-2556 


27.  Kondratyev,  K.  Ya.  and  Yu.  M.  Timofeev,  "Numerical  Ex- 
periments on  Studying  Transmission  Functions  of  Atmos- 
pheric Gases,"  in  Radiation  Including  Satellite  Tech- 
niques.  Proceedings  of  WMO/IUGG  Sumposium,  WMO  Tech. 

Note  104  (1968)  . 

28.  Armstrong,  B.  , "Analysis  of  tne  Curtis-Godson  Approxi- 
mation and  Radiation  Transmission  Through  Inhomogeneous 
Atmospheres,"  J.  Atmos.  Sci. , Vol.  25  (1968),  p.  312. 

29.  Zdunkowslci,  W.  and  W.  Raymond,  "Exact  and  Approximate 
Transmissio.  Calculations  for  Homogeneous  and  Nor.- 
Homogeneous  Atmospheres,"  Beitr  zur  Phys.  der  At  nos., 

Vol.  43  (1970)  , p.  185. 

30.  Cantor,  D.  and  J.  Evans,  "On  Approximation  by  Positive 
Sums  of  Powers,"  SIAM  J.  Appl.  Math,  (1970),  Vol.  18, 
p.  380. 

31.  PenjtCorf , R.  , "Tables  of  the  Refractive  Index  for  Stand- 

ard Air  and  the  Rayleigh  Scattering  Coefficient  for  the 
Spectral  Region  Between  .2  and  20  d.  Opt.  Soc.,  Amer. , 

Vol.  47  (1957),  p.  176. 

32.  Peck,  E.  and  K.  Reeder,  "Dispersion  of  Air,"  J.  Opt.  Soc.  , 
Amer. , Vol.  62  (1972) , p.  958. 

33.  Gucker,  F.  , et.  al.,  "Intensity  and  Polarization  of  Light 
Scattered  by  Some  Permanent  Gases  and  Vapors,"  J.  Chem. 
Phys. , Vol.  50  (1969),  p.  2526. 

34.  Dave,  J.  , "Subroutines  for  Computing  the  Parameters  of 
the  Electromagnetic  Radiation  Scattered  by  a Sphere," 

IBM  Program  360D-17. 4002 , IBM  Scientific  Center,  Palo 
Alto,  Calif.  (1968). 

35.  Kattawar,  G.  and  G.  Plass,  "Flectromagnet ic  Scattering 
From  Absorbing  Spheres,"  Appl.  Opt.,  Vol.  6 (1967), 

p.  1377. 

36.  Denman,  H.  , W.  Heller,  and  W.  Pangonis,  Angular  Scatter- 
ing Functions  for  Spheres,  Wayne  State  Univ.  Press,  De- 
troit  (1966)  . 

37.  Deirmendj ian , "Tables  of  Mie  Scattering  Cross  Sections 
and  Amplitudes,"  Rand  Report  R-407-PR,  Rand  Corp. , Santa 
Monica,  Calif.  (1963). 


R-3 


SSS-R-75-2556 


38. 


39. 


40. 


41. 


42. 


43. 


44. 


45. 


46. 


47. 


48. 


49. 


50. 


51. 


Irvine,  W.  and  J.  Pollack,  "Infrared  Optical  Properties  of 
Water  and  Ice  Spheres,"  Icarus,  Vol.  8 (1968),  p.  324. 


Rusk,  A.,  D.  Williams,  and  M.  Querry,  "Optical 
of  Water  in  the  Infrared,"  J.  Opt.  Soc. . Amer . . 
(1971),  p.  895.  


Constants 
Vol.  61 


Robertson,  C.  and  D.  Williams,  "Lambert  Absorption  Coeffi- 

SOC.,  A.er., 


Querry,  M. , et.  al.,  "Optical  Constants  in  the  Infrared 
for  Aqueous  Solutions  of  NaCl,"  o.  Opt.  Soc.,  Amer., 
Vol.  62  (1972),  p.  349.  ^ ' 


Twitty,  J.  and  J.  Weinman,  "Radiative  Properties  of  Car- 
bonaceous Aerosols,  " J.  Appl.  Met..  Vol.  10  (1970),  p.  V25. 


Dalzell,  W.  and  A.  Sarofim,  "Optical  Constants  of  Soct  and 
Their  Application  to  Heat-Flux  Calculations,"  J.  Heat  Trans- 
fer, (February  1969),  p.  100.  


Spitzer , 
Quartz, " 


W.  and  D.  Kleinman,  "Infrared  Lattice  Bands  of 
Phys.  Rev.,  Vol.  121  (1961),  p.  1324. 


Volz,  F.,  "Infrared  Absorption  by  Atmospheric  Aerosol  Sub- 
stances," JGR,  Vo.  77  (1972),  p.  1017. 


Volz,  F.,  "Infrared  Refractive  Index  of 
Substances,  •' Appl.  Opt.,  Vol.  11  (1972), 


Atmospheric  Aerosol 
p.  755. 


Blanco,  A.  and  G.  Hoidale,  "Infrared  Absorption  Spectra  of 
Atmospheric  Dust,"  ECOM-5193,  1968,  Atmospheric  Sciences 
Lau,  Lnite  Sands  Missile  Range,  New  Mexico. 


Flanigan,  D.  and  H.  DeLong,  "Spectral  Absorption  Character- 
istics of  the  Major  Components  of  Dust  Clouds,"  EATR-4430 
1970,  Edgewood  Arsenal,  Maryland.  ' 


Hunt,  J.,  M.  Wisherd,  and  L.  Bonham,  "Infrared  Absorption 
Spectra  of  Minerals  and  Other  Inorganic  Compounds,"  Anal. 
Chem.  , Vol.  22  (1950),  p.  1478. 


Went,  F.,  "Cn  the  Nature  of  AitJeen  Condensation  Nuclei," 
Tellus,  XVIII  (1966),  p.  549. 


Junge,  C.  E. , Air  Chemistry  and  Radioactivity,  Academic 
Press,  New  York  (1963)  . ^ 


R-4 


SSS-R-75-25b6 


52.  Deirmendj ian,  D.  , Electromagnetic  Scattering  on  Spherical 
Polydispersions , Elsevier,  New  York  (1969)  . 

53.  Khrgian,  ed. , Cloud  Physics,  Israel  Program  for  Scientific 
Translations,  Jerusalem  (1963) , p.  83. 

54.  Germogenova,  0.  ^.,  "Atmospheric  Haze:  A Review,"  Report 

1821,  Bolt  Deranek  and  Newman,  Cambridge,  Mass.  (1970). 

55.  Twomey,  S.,  H.  Jacobowitz,  and  H.  Howell,  "Light  Scattering 
by  Cloud  Layers,  J.  Atmos.  Sci.,  Vol.  24  (1967)  , p.  70. 

56.  Blifford,  I.  H. , "Tropospheric  Aerosols,"  JGR,  Vol.  75  (1970), 
p.  3099. 

57.  Chang,  D.  and  R.  Wexler,  "Meteorological  Relationships  to 
Tropospheric  and  Stratospheric  Turbidity  Profiles,"  AFCRL- 
70-0197,  Air  Force  Cambridge  Research  Labs.,  Cambridge, 

Mass . (1970)  . 

58.  Kondratyev,  K.  Ya.,  et.  al.,  "Aerosol  Structure  of  the  Tro- 
posphere and  Stratosphere,"  WMO  Technical  Note  No.  104, 
Radiation  Including  Satellite  Techniques  (1968) . 

59.  Dufour,  L. , "The  Atmospheric  Aerosol,"  translated  at  Atmos- 
pheric Sciences  Laboratory,  White  Sands  Missile  Range,  Now 
Mexico  (1969)  . 

60.  Selezneva,  E.  S.,  Atmospheric  Aerosols,  translated  and  pub- 
lished for  ESSA  (now  NOAA)  by  Indian  National  Scientific 
Documentation  Center,  New  Delhi,  India  (1970) . 

61.  Elte  .'mann,  L.  , "Parameters  for  Attenuation  in  the  Atmos- 
pheric Windows  for  Fifteen  Wavelengths,"  Appl.  Opt.,  Vol.  3 
(1964) , p.  745. 

62.  Meszaros,  A.,  "On  the  Variation  of  the  Size  Distribution 
of  Large  and  Giant  Atmospheric  Particles  as  a Function  of 
the  Relative  Hamidity,"  Tellus , YXIII  (1971),  p.  436. 

63.  Hanel,  G. , "The  Real  Part  of  the  Mean  Complex  Refractive 
Index  and  the  Mean  Density  of  Samples  of  Atmospheric  Aero- 
sol Particles,"  Tellus , XX  (1968),  p.  371. 

64.  Ralston,  A.  and  H.  Wilf,  eds..  Mathematical  Methods  for 
Digital  Computer s , Vol.  II  (1967) , Wiley,  New  York. 

65.  Dave,  J. , "Effect  of  Coarseness  of  the  Integration  Incre- 
ment on  the  Calculation  of  the  Radiation  Scattered  by  Poly- 
dispersed  Aerosols,"  Appl . Opt. , Vol.  8 (1969),  p.  1161. 


R-5 


SSS-R-75-2556 


66. 


67. 


68. 


69. 


70. 


71. 


72. 


73. 


74. 


75. 


76. 


77. 


78. 


79. 


Dave,  J.,  "Effect  of  Varying  Integration  Increment  on  the 
Computed  Polarization  Characteristics  of  the  Radiation 
Scattered  by  Polydisper sed  Aerosols,"  Appl.  Opt.  , VoJ  . o 
(1969)  , p.  2153. 

Hansen,  J. , J.  Atm.  Sci . , Vol.  26  (1969),  p.  478. 

Potter,  J.  , "The  Delta  Function  Approximation  in  Radiative 
Transfer  Theory,"  J.  Atm.  Sci. , Vnl.  27  (1970),  p.  943. 

Grant,  I.  and  G.  Hunt,  "Solution  of  Radiative  Transfer 
Problems  Using  the  Invariant  S^  Method,"  Mon.  Not.  Roy^ 

Astron.  Soc . , Vol.  141  (1968),  p.  27. 

Grant,  I.  and  G.  Hunt,  "Discrete  Space  Theory  of  Radiative 
Transfer  I.  Fundamentals,"  Proc.  Roy.  Soc.  Lond.  A.,  Vol. 

313  (1969) , p.  183. 

Grant,  I.  and  G.  Hunt,  "Discrete  Space  Theory  of  Radiative 
Transfer  II.  Stability  and  Non-Negativity,"  Proc.  Roy.  Soc. 
Lond.  A. / Vol.  313  (1969),  p.  199. 

Grant,  I.  and  G.  Hunt,  "Solution  of  Radiative  Transfer  Prob- 
lems in  Planetary  Atmospheres,"  Icarus,  Vol.  9 (1968),  p.  5Zb. 

Hunt,  G.,  "The  Effect  of  Coarse  Angular  Discretization  on 
Calculations  of  Radiation  Emerging  from  A Model  Cloudy  At- 
mosphere," JQSRT,  Vol.  11  (1971),  p.  309. 

McCleese,  D. , J.  Margolis,  and  G.  Hunt,  "Laboratory  Simula- 
tion of  Absorption  Spectra  in  Cloudy  Atmospheres,  Nature 
(Phys . Sci . j , Vol.  233  (1971),  p.  102. 

Preisendorfer,  R.  W. , Radiative  Transfer  on  Discrete  Spaces, 
Pergamon  Press,  Oxford  (1965) • 

van  de  Hulst,  H.  C.  , "A  New  Loolc  at  Multiple  ^^catterlng , " 
unpublished  repcit  (1963),  NASA  Inst,  for  Space  Studies, 

New  Yorlc. 

Rybiclci,  G.  and  P.  Usher,  "The  Generalized  Ricatti  Trans-^^ 
formation  as  a Simple  Alternative  to  I:i/ariant  Imbedding, 
Astrophys.  J. , Vol.  14C  (1966),  p.  871. 

Lathrop,  K.  and  B.  Carlson,  "Numerical  Solution  of  the 
Boltzmann  Transport  Equation,"  J.  Comp.  Phys.,  Vol.  2 (1967), 
p.  173. 

Kroolc,  M.  "On  the  Solution  of  Equations  of  Transfer  I," 
Astrophys.  J.,  Vol.  122  (1955),  p.  488. 


R-6 


SSS-R-75-2556 


80.  Gates,  W.  L.  , et.  al.,  "A  Documentation  of  the  Mintz- 
Arakawa  Two-Level  Atmospheric  General  Circulation  Model," 
Report  R-877-ARPA,  Rand  (December  1971). 

81.  Green,  A.,  Appl.  Opt. , Vol.  3 (1964),  p.  203. 

82.  Weller,  G.,  S.  Bowling,  K.  Jayaweera,  T.  Ohtaki,  S.  Parker, 
G.  Shaw,  and  G.  Wendler,  "Studies  of  the  Solar  and  Terres- 
trial Radiation  Fluxes  Over  Arctic  Pack  Ice,"  (1972),  Geo- 
physical Institute,  University  of  Alaska,  Fairbanks,  Alaska. 

83.  Valley,  S.,  ed.  , Handbook  of  Geophysics  and  Space  Environ- 
ments, McGraw-Hill,  New  York  (1965). 

84.  Huschke,  R. , "Arctic  Cloud  Statistics  from  'Air-Calibrated' 
Surface  Weather  Observations,"  RAND  Memorandum  RM-6173-PR 
(1969),  The  RAND  Corporation,  Santa  Monica,  California. 

85.  "Inadvertent  Climate  Modification;  Report  of  the  Study  of 
Man's  Impact  on  Climate  (SMIC),"  (1971),  MIT  Press,  Cam- 
bridge, Mass. 


R-7 


